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FITTING MULTIDIMENSIONAL SPLINES USING STATISTICAL VARIABLE 

SELECTION TECHNIQUES 

By 

Patricia L. Smith ^ 

SUMMARY 

This report demonstrates the successful application of statistical 
variable selection techniques to fit splines. Major emphasis is given to 
knot selection, but order determination is also discussed. Tvo FORTRAN 
backward elimination programs using the B-spline basis were developed, and 
the one for knot el aination is compared in detail with two other spline- 
fitting methods and several statistical software packages. An example is 
also given for the two-variable case using a tensor product basis, with a 
theoretical discussion of the difficulties of their use. 

1. INTRODUCTION 


Polynomial splines have often been employed in modeling or data fitting 
when the functional form of the relationship between the dependent and inde- 
pendent variables is unknown. The major problem has been how to avoid 
under- or overfitting the data. A strictly mathematical approach is to add 
knots one at a time and move them around until the L 2 (or some other) norm 
of the errors is less than a preselected tolerance level (ref. 1). A major 
problem with this approach is that a good fit depends entirely on the suo- 
jective selection of the tolerance level. A fitting method which attempts 
to avoid this problem is the smoothing technique introduced by Reinsch 
(ref. 2), but it requires the experimenter to have good a priori information 
■ibout the data 0 - Che process which generated it. Both of these methods are 
currently feasible only for functions of a single variable. 

A statistical approach to the curve-fitting problem using the method of 
cross-validation was introduced by Wahba and VIold (ref. 3). The major 

^Assistant Professor, Department of Mathematics, Old Dominion University, 
Norfolk, Virginia 23508. 



advantage of this procedure is its automation: no a priori information is 

needed. There are several disadvantages, however. Every data point is a 
knot so that the resulting functional form is difficult to use and inter- 
pret. In addition, if there are clearly identifiable trends in certain 
portions of the data such as linearity or sharp bends, this information is 
lost analytically even though it show up when the spline is plotted. The 
practical use of this technique is also currently restricted to functions of 
one or two variables. The two variable case is considered in Wahba (ref. 

4), with higher dimensions discussed in Wahba and Wendelberger (ref. 5). 

Other statistical approaches to the variable knot spline problem have 
considered the knots as parameters in the model. However, this presents 
problems in finding the least squares solution and in subsequent statistical 
estimation and testing procedures because the model is nonlinear. Tradi- 
tional (ref. 6) as well as Bayesian (ref. 7) approaches have been investi- 
gated, but both are limited in scope and application. Further, in most 
cases, though the knot locations have been variable, their number has been 
fixed £ priori by the analyst. Some exceptions are the works of Ertel and 
Fowlkes (ref. 8), Smith and Smith (ref. 9), and Agarwal and Studden (ref. 

10), but, as with most other approaches mentioned above, they have not been 
developed to fit splines in several variables. 

The technique investigated in this research is the use of variable 
selection procedures to fit splines. If a pool of knots is fixed in advance, 
then statistical linear models theory can be applied in a variable selec- 
tion framework. There are four major advantages of the variable selection 
approach to fitting splines. First, variable selection procedures are 
essentially user independent (automatic) in their use of the F test as a 
stopping criterion. Second, they are widely available in statistical soft- 
ware. Third, final fit’s may have straightforward intepretations because of 
their simplicity or theoretical foundation. Fourth, regression diagnostics, 
such as outlier detection, may be performed. These advantages and other 
desirable properties are discussed in Section 4, along with a comparison of 
several methods and software. 

The theory applies not only to splines in a single variable, but also 
to splines in several variable using a tensor product basis. However, as 
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the careful and detailed development of this technique in the one-variable 
case is considered a crucial step to its use in several variables, discus- 
sion of the multivariate case is restricted to Section 7, and includes an 
example of its successful application to aerodynamic modeling. 

The major emphasis of this report is the application of variable selec- 
tion procedures for choosing the number and locacioii of the knots for 
splines in a single variable of fixed order (“ degree +1). A detailed dis- 
cussion of this "knot selection" approach is given in Section 2 with exam- 
ples, comparison of methods and software, and applications in Sections 3-5. 
Choosing the spline order with the number and location of the knots fixed is 
of less interest and considered in Section 6 only. FORTRAN programs which 
apply backward elimination in these two contexts were written as part of 
this research and discussed in Sections 2 and 6. Their documentation, flow- 
charts, and listings are given in the Appendix. 
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2. THE KNOT SELECTION (KS) PROCEDURE 


Statistical variable selection procedures can be used as a KS procedure 
to chcoKe the number and location of knots in fitting splines. The "+" 
function basis is suitable for this, at least theoretically, because it is 
easily interpreted. Knots and knot multiplicities correspond to individual 
terms so that selection or deletion of terms is equivalent to sel ction or 
deletion of knots. The knots are thus selected indirectly. For example, a 
continuous linear spline with knots t ,...,t may be written as 
I 

B •*■6, x+E-B.Cx-t.)^, where u , * u for u > 0 '•nd zero other- 
o 1 2 1 1 +’ + 
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wise. 


Selection of the "spline term" (x - i® actually selection of the 

knot t^. Because we don't know where the breakpoints should be, we provide 
as candidate variables a liberal number of spline terms, i.e,, a pool of 
knots, more > “ expect or want to eventually use, and blanket the 

domain. iS , the ’tual number and location of the knots used in the final 
model is LMiKnown at the beginning in the sense that we are selecting from a 
larger set. 

While "+" functions are easily defined in current statistical software 
packages and fit into the statistical hypothesis testing framework without 
modification (ref. 11), computational problems such as carry-over in round- 
off error and mult icoll inear ity greatly restrict their use. As will be seen 
in Section 4, the backward elimination (stepdown) procedures are especially 
troublesome because all terms must be fit initially. An alternative is the 
use of the computationally advantageous B-spline basis (ref. 1). Unfortu- 
nately, it does not fit ea.ily into the hypothesis testing framework and 
cannot be used in existing statistical software packages. There was thus a 
need for the development of a KS procedure using B-splines. Construction 
of hypotheses which are useful in B-spline regression, including testing the 
importance of knots, has been detailed in Smith (ref. 12). As part of this 
research, these results have been implemented in two FORTRAN computer 
programs, one of which accommodates the backward elimination of knots using 
the B-spline basis. Examples in Section 3 give the results of using this 
FORTRAN program, and comparisons with several statistical software packages, 
as well as with other statistical spline-fitting methods, are detailed in 
Section 4. 

The use of variable selection is a sort of compromise between the tech- 
niques which use either fixed or variable knots. Its most important advan- 
tage, and one which makes possible all others, is that because the maximum 
number and location of the knots is fixed in advance, the statistical theory 
of general linear models applies. Consequently, the least squares solution 
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is easily obtained at any given step, and hypothesis testing and interval 
estimation are straightforward. As mentioned earlier, details for using the 
B~spline basis are given in reference 12. The selection of knots can thus 
be accr/mplished through t tests. This fits exactly into the variable 
selection framework for (1) spline models in a single variable, (2) models 
in several variables with spline terms in one or more variables, and (3) 
models in several variables with tensor products defining higher dimensional 
splines. Also, trends in the data in one or more variables may be easily 
detected through the selection of a few knots. Several examples of this 
will be given in the next section. Further, in some experimental situa- 
tions, models may be easily interpreted because the coefficients are physi- 
cally meaningful, as in some examples in Sections 5 and 7. 

3. EXAMPLES (ff THE KS PROCEDURE 

Four data sets were examined using the FORTRAN knot elimination pro- 
gram. The maximum number of continuity constraints allowed for any given 
order were imposed. The first data set, the Indy data, is rather simplistic 
but has appeared in the statistical literature several times in connection 
with curve-fitt > i splines. It is a record of the average winning 
speeds at the L d .polis 500 from 1911-1971, except for 1917-1918 and 
1942-1945, during the two World Wars when the race was not run. Poirier 
(ref. 13) fit the data with a cubic spline with 2 knots, one each at the 
midpoint of the non-racing years. The data were coded so that x ■ year - 
1910 with knots 7.5 and 33.5. The output and graphs from the knot elimi- 
nation routine are shown in Figures 3.1 to 3.4, with circles around the 
function values of the l:nots. Using an F-table value of 8.0 (a ■ 0.01), the 
KS procedure eliminates both knots so that a cubic polynomial is adequate 
to fit the data. If a linear rather than a cubic spline is fit, only the 
knot at X ” 7.5 can be eliminated (Figures 3.5 to 3.7). 

The second example is noisy data generated from the function used in 
reference 3 

f(x) = 4.2b(e - 4e + 3e ) 
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Figure 3.1. Ou' p^it for knot elimination. 


Indy data. Cubic spline. 
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for jc[0,3]. For x starting at zero, we generated 100 data points at 
intervals of 1/32 up to 99/32 and added normal random noise, N(u “ 0, 
a * .2), the value of o the same as that used by Wahba and Wold (WW) . A 
graph of the function and generated data is shov in Figure 3.8. Figures 
3.9 to 3.27 show graphical results of the si’epc /n procedure for cubic 
splines starting with 19 equally spaced interi?i knots, and using an F-table 
value of 8.0. By examining this sequence of g.aphs, it becomes clear how 
the elimination of knots makes the spline smoother by making it less noise 
dependent . 

An F-table value of 4.0 (a >« 0.05) rather .han 8.0 results in stepdown 
terminating with 5 knots remaining (Fig. 3.23, >. 20). The latter fit is 
more data dependent and clearly inferioi i.n tci\as of recovering the desired 
function. Use of the larger F value thus seems appropriate and keeps the 
procedure from terminating "prematurely.” Graphs of starting and ending 
fits to the data, beginning with 39 interior knots, are shown in Figures 
3.28 to 3.29, and the results are roughly the sane as when 19 knots are used 
initially (Figure 3.27, p. 22). A phenomenon which occurs throughout most 
of these fits is the downward hook in the upper range of the x's due to a 
cluster of 3 data points. Figure 3.30 shows the conclusion of stepdown with 
those 3 points omitted and helps to illustra .e the fact that different noise 
results in different fits. 

Tlie method used by Wahba and Wold to recover the function is a modifi- 
cation of the smoothing technique introduced by Re insch (ref. 2). They use 
cross-validation to determine the smoothing paraneter, and their resulting 
fit is shown in Figure 3.31. Referring again to Figure 3.27, r. . 22, we see 
that the results of the two methods compare very f ivorably A more detailed 
comparison of these methods and others is made in the r ' vt section. 

Smith and Smith (SS) (ref. 9) examine a scale^ version of the WW 
function, f(x) * 4.26 (e - 4e + 3e for xe[o,l]. A sample 

of size 600 equally spaced points was generfed, ind a variance of 0.039 (as 
in SS ) was used for the normally distrib- ted zero mean noist. Results from 


12 



HU DflTP 


Figure 3.8 


Figure 3. 


. The Wahha-Wold (WK) function and data generated from it 



u. First stop of knot elimination. W3'/ data. 


WW DflTfi 


>8 


3. 


Figure 


'^igure 


.10. Second step of knot elimination. WW data. 



.11. Third step of knot elimination. W data. 
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Figure 3.12. Fourth step of knot elimination. WW data. 
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Figure 3.13. Fifth step of knot elimination. WW data. 
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Figure 3.14, Sixth step of knot elimination. WW data. 



Figure 3.15. Seventh step of knot elimination. WW data. 
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‘ Figure 3.16. Eighth step of knot elimination. data. 



Figure 3.17. Ninth step of knot elimination. WW data. 
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Figure 3.18. Tenth step of knot elimination, WU data. 



Figure 3.19. Eleventh step of knot elimination. WW data 
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Figure 3.20. Twelfth step of knot elimination. WW data. 
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Figure 3 21. Thirteenth step of knot elimination. WW data. 
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Figure 3.22. Fourteenth, step of knot elimination. WW data. 



Figure 3.23. Fifteenth step of knot elimination. WW data. 
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Figure 3.30. Final step of knot elimination with 3 data points in 
upper x-range omitted- WW data. 
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Figure 3.31. Spline fit obtained by cross-validation by Wahba and 
Wold. WW data. 



twit stopdi^i%m runs fitting cubic splines are ahowti in Figurea 3.32 to 3.33, 
beginning with l<) and 44 knots. Three knots i-enain in Figure 3.32 with a 
slightly wigglier fit than that in Figure 3.33 with one remaining knot. 

Those results show that a larger knot selection pool allows reduction to 
possibly a fewer number of final knots and a smoother fit, which, for 
simplicity, is more desirable. 

Smith and Smith use asymptot ic results to determine a stoppit\g rule for 
adding knots one at a time tvi the mvxlel. Figure 3.34 shv^ws their results 
using cubic splines overlaid on the true functiv'n. The data were not plot- 
ted svi that the distinctions betiwen the two functions would not be lost. 
Applying stepiiown using these 4 initial knv>ts resulteil in Figure 3.3S, a fit 
which mnooths the wiggles visible in Figure 3.34. As seen in the two 
previous f igvires , hv'wever, using a larger pool vif knots results in a 
»u 20 other and more satisfactory recovery of the function. The SS method is 
cvxnparevl in more detail to both the WV) and KS meth^xls in the next 
sec c ion. 

The final function examinevl is f(x) * sin (x*3 for xc[o,4.5], which 
allk'ws tor more than tw.i periods of the sine w.sve and gradually increases 
the frequency. Three hundred data points were used with o • ,2 for the 
norra.s! noise. Beginning and ending cubic spline fits frv>m a stepdown run 
are shown in Figures 3.3b to 3.37, starting with 14 intei 1. ; knots and 
ending with 4 . We note that more knots .are needed for the final fit than 
foi the functions previously discussed due to the increased curvature of the 
function. Moat of the wiggliness in the initial spline fit occurs on the 
more gradual slope at the lower end of the x-range and is removed as knots 
are removed. This phenomenon also occurs on the "flat" pi^rtion of the SS 
and UW data. 

In order to .assess the effects of a lower noise level on the KS tech- 
nique, randiMxi variables used for the noise on the WW function were generated 
using 0 " .1 anil .05. Final fits are shv'wn in Figures 3.38 to 3.34, and 
referring back to Figure 3,27, p. 22, which shows results using o ■ .2, we 
see that titring data with a lower noise level results in uKare knots 
remaining at the end of the prc>cedure. This tendency is especially striking 
when d.ita fr.w the function itselt is fit, that is, when no noise is added 
so that to 
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.36. First step of knot elimination with 19 knots, 
function is sin (x^). 



.37. Final step of knot elimination from 19 knots, 
function is sin (x^). 
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True 
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recover the function we actually need to interpolate. A stepdown from 19 
knots results in a spline with 12 knots as shown in Figures 3.40 to 3.41. 
Both the true and fitted functions are graphed, but there is no perceptible 
difference between the two. 

Wiggliness in data should be smoothed (i.e., ignored) if it is per- 
ceived as noise, but should be fit if it is perceived as trends in the 
underlying process. Thus, a danger in applying the KS technique is using 
too small or too large a pool o.' knots, the former problem is illustrated 
quite well in Figures 3.42 to 3. +3, where noisy data generated from sin (x^) 
is fit with the KS technique beginning with too few knots to allow the 
bending necessary to recover the function, especially near the third peak. 

It is interesting to see that the three knots eliminated were in the lower 
end of the x range where the underlying function is not wiggly. A 
comparison of Figures 3.37, p. 28, and 3.42 reveals that both have 9 knots, 
but a better fit is obtained from the one which began with 19 knots (Fig. 
3.37): its 9 knots are more selectively and better placed. 

4. COMPARISON OF METHODS AND SOFTWARE 


In the previous section, two functions introduced in the literature (WW 
and SS) were examined using the FORTRAN knot elimination program. The pur- 
pose was to compare results, which we do in this section, in light of what 
we consider to be the most desirable properties of curve-fitting with 
splines, these are: 

( 1) good results ; 

(2) computational efficiency; 

(3) diagnostics capabilities; 

(4) user independence; 

(5) ease of interpretation; and 

(6) ease of use. 

We also give in this section the results of using several statistical soft- 
ware packages on the Indy and WW data, fitting both linear and cubic 
splines. 
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Most statisticians have ready access to variable selection procedures, 
either in programs they have written themselves, or in widely available 
statistical software packages. Fitting splines through knot selection with 
these programs is a potential advantage of their use, which is realized only 
if good results are obtained. A summary of the results of using four such 
packages is given in Table 4.1: SAS (ref. 14), SPSS (ref. 13), MINITAB 

(ref. 16), and BMDP (ref. 17). 

Table 4.1. Results of using variable selection techniques to fit 
splines with four statistical software packages. 






Stepwise 


Stepdown 





Indy 

WW^ 

Indy 

WW^ 

SAS 

linear 



✓ 

/ 

✓ 

✓ 


cub ic 



✓ 

/ 

/ 


SPSS 

linear 



/ 

✓ 

/ 

✓ 


cub ic 



/ 

✓ 


✓ -2 

MINITAB 

linear 



/ 

✓ 

✓ 

✓ 


CL'b ic 



✓ 

X 

X 

X 

BMDP^ 

linear 



X 

X 

X 

X 


cub ic 



X 

X 

X 

X 

^Selection pool of 

19 

interior 

knots . 




^Numerical output 

has 

some inaccuracies, but 

overall 

results are 

correct. 

^Tolerance cannot 

be 

made low 

enough to force 

entry 

of necessary 

terms. 


In the case of stepwise procedures, accuracy was determined oy comparing 
outputs for the various packages among themselves, while outputs for the 
stepdown procedures were compared with the FORTRAN B-spl ine knot elimination 
program. Results are surprisingly good considering the fact that the "+" 
function basis must be used. Entries marked with an "X" indicate failure to 
produce accurate results or, sometimes, any results at all due to high 
mul ticoll inearity in the models or low tolerance, especially in stepdown. 
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The minimum tolerance allowed for SPSS, 10 “^ 2 ^ j,g ^g force entry 

of some of the polynomial terma or Co get results in stepdown. For BMDP, 
the tolerance of O.Ol for variable selection was not low enough to force 
entry of necess y terms Co get results for any of the cases considered. As 
expected, less trouble was had with fewer knots (Indy data), lower .egree 
(linear), and simpler models (stepwise^. Stepdown gave accurate results in 
several cases even for a large niinber of knots, but there are limitations. 
For instance, computational problems were encountered by SAS for the cubic 
WW data with 39 knots. The final odels determined by s*'epwise and step- 
down, however, were either identical or very similar. The occasional user 
of splines could thus safely rely on stepwise procedures from one of several 
packages to give good results. 

Table 4.2 compares several spline-fitting methods: Wahba-Wold (WW), 

Smith-Smith (SS), and knot selection (KS). As the latter method may be 
implemented through several different computer progra-iis, two statistical 
packages and the FORTRAN knot elimination routine are included. All methods 
give good results for the data examined, though as seen in earlier discus- 
sion, care must be taken when using the statistical packages, especially 
for stepdown. Their use of the "+" function makes them computationally 
inefficient and can cause severe problems. They are handy, however, for the 
occasional user as is the WW method which is available as an IMSL subroutine 
(ref. 18). The KS techniques depend on setting an a level for the hypo- 
thesis tests and specifying an initial pool of knots but are otherwise user 
independent. The WW method is "completely automatic," while the SS method 
depends on user application of the stopping criterion. The KS approach in 
general produces results which are >ier to interpret. 

Results from this section and from Section 3 show that splines 
fit by knot selection recover the underlying functions quite well and 
compare very favorably with the results of Wahba and Wold and improve upon 
those of Smith and Smith. Though somewhat simplistic, the knot selection 
approach provides an alternative to the method of cross-validation and 
offers a great computational savings. In addition, there is the possibility 
of analytic or physical interpretation in many modeling situations, an 
exanple of which is given in the next section. 
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Table 4.2. Comparison of 
and software. 

spl ine- 

fitting techniques 






Knot Selection 

Desirable Properties 



SAS 

SPSS 

B-Spl ines 
FORTRAN 

Good results 

/ 

/ 

/ 

/- 

✓ 

Computational efficiency 

✓ 

✓ 

X 

X 

✓ 

Diagnostics capabilities 

X 

X 

/ 

✓ 

✓ 

User independence 

✓ 


✓ - 

✓ - 

✓- 

Ease of interpretation 

X 

X 

✓ 

/ 

✓ 

Ease of use occasionally 

✓ 

X 

/-I 

/ 

X 


^ SAS is available only on IBM-compatible machines. 


3. SOME SPECIAL APPLICATIONS 

Probably the most useful .•’ppl icat ion of the KS technique is data- 
smoothing, and in Section 3 we saw several examples of recovering underlying 
functions from noisy data. A v:>riacion that is useful in simulation 
experiments is smoothing the ample quantile function (ref. 19). This 
function is a left-continuous step function defined as Q(u) ■ for 

(i-l)/n < u < i/n, where n is the sample size and fs the i-th 

order statistic. Experimental conditions can be simulated by generating 
data which behaves like the original, and a smoothed sample quantile 
function provides a continuous distribution from which to draw the simulated 
data. An advantage of smoothing the sample quantile function, rather than 
its pseudo-inverse, the sample cumulative distribution function, is that the 
former always has domain [o,l] regardless of the type of distribution. 
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Progiamming can thus b> standardized, as, for example, in the determination 
of the original knot selection pool. 

The KS technique is also useful in modeling. For example, stepwise 
regression has been applied successfully by Klein, Batteraon, and Smith 
(ref. 20) to model flight data using splines. They use function terms 
defined in the angle-of-attack variable in a Taylor series expansion of 
force and raoraeni; coefficients in order to model longitudinal motion of an 
airplane. CXie of their simple "spline-modified" Taylor series expansiono of 
the vprtical aerodynamic force coefficient ig given by 


C 

7 . 


C (a) , + C (a) 

Z q "0 z 


Z q 
6 “0 


+ C (a) 6 
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where 


C (a) - C (a 
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C (a) 
^6 


- C 
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+■ 


U3 

I 

1-2 


a - Sjj ) 


0 

+ 


and a is the angle o^ attack, q' is the nondiitensional pitch rate, 6^ 

is the elevator deflection, C ■ 3C /3a, C *■ 3C /3q', C ■ 3C /3^ . 

z z z z z> z e 

a q 0 

^ e 

They then use stepwise regression to select terras, and thus knots, in the 
model. Thi;j spline representation preserves the concept of stability and 
control derivatives inherent in the usual Taylor series expansion of ae’"o- 
dynamic coefficients but has the advantage of providing a representation of 
C over an extended range of the angle of attack a. A global model over 
the observed range of a is thus obtained through the use of splines. 
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6. OTHER USES OF VARUBLE SELECTION PROCEDURES IN SPLINE REGRESSION 


Thus far we have emphasized the use of variable selection to choose the 
nvanber and location of knots. There are other possible, but perhaps less 
useful, "extensions" to spline regression of variable selection procedures 
based on polynomial or multiple regression models. In the latter cases, the 
purpose is to d termine the polynomial degree and the important independent 
variables anc .nteractions. This is accomplished by examining the 
contribution of individual terms in the model. With univariate spline 
models, however, there are several polynomial pieces, not just one, whose 
degrees may be examined, and, as seen previously, we may examine the 
importance of each knot. Also, one may wish to examine the continuity 
conditions at one or more breakpoints as in the example discussed by Smith 
(ref. ll). Thus, the complexity of the spline model over the polynomial 
model manifests itself in the greater number of ways the dimension of the 
spline parameter space may be altered. %>lines in several variables present 
even more possible diversity since, for example, t*K>-var iable spline 
continuity occurs not across points but along lines connecting grid 
po ints. 

While it might be nice to have a single software package which could 
perform any combination of these spline hypothesis tests, it is neither 
feasible nor desirable. The major i ason is that variable order splines, 
i.e., splines with polynomial pieces of different degrees, have not been 
sufficiently researched by mathematicians to allow for the satisfactory 
construction in a general framework of a basis using either functions or 
B-splines. Lowering or raising the degree of a single polynomial piece must 
be accomplished by applying restrictions to the model, and hypothesis tests 
must then use restricted least squares. In simple cases this may be 
straightforward (references ll and 21), but in general the task is unmanage- 
able. For example, the u r is subject to hidden analytical errors as when 
the regression or hypothesis degrees of treedom are not equa^ to the number 
of restrictions because some restrictions are obtained automatically through 
linear combinations of others. While theoretically such dependencies can be 
checked, the usual methods would need some revision in the case 



of B-spline regres ion since hypotheses involve values of the fitted spline 
or its derivatives (ref. 12). In the case of the function basis, most, 
but not all, of the individual terms are meaningful. Hoifever, the innocent 
yet indiscriminant selection or removal of terms through hypothesis tests 
can result in fits which are statistically valid yet nonsensical because 
they are uninterpretable in terms of polynomial degree or knot locations 
(ref. 11). Because of these various difficulties, it is reasonable to con- 
struct task-specific procedures. 

The application of variable selection to knot selection, as in the 
examples in Section 3, is useful for smoothing data with a fixed order 
spline with maximian continuity conditions. In these cases the interest is 
not in the spline order but rather in determining the minimal number of 
knots deemed adequate to faithfully represent the data. Cubic splines are 
popular because of their low degree and second derivative continuity. The 
selective use of forward or backward algorithms in some statistical software 
packages using functions (see Section 4), or the backward elimination 
FORTRAK program developed here using B-splines, nay be used for this 
purpose. 

Another possible "extension" of variable selection to splines is the 
determination of the polynomial degree while keeping the number and location 
of knots fixed, that is, not consider the knots as "variables" to be either 
entered or removed. Because of the difficulties with variable order splines 
discussed above, we must restrict ourselves to polynomial pieces of the sane 
degree. Unfortunately, even further constraints are necessary for this 
version. The ideal situation would be to compare a maximally continuous 

[c )k-th order spline, i.e., a k-th order spline with continuous f, 

(1) (k-2) ... . r 

t ,...f , with a maximally continuous k-l-st order spline (C J. 

A formal test, however, is not possible. This can be easily seen by consid- 
ering a specific example using the partial ordering of some spline models 
given in reference 11. Basis elements for C® and quadratic splines 
and for CP linear splines with one knot are shown in Fig. 6.1. A compari- 
son of orders 3 and 2 (degrees 2 and 1) which retained maximum continuity 
conditions would require comparing the quadratic with the C® linear. 
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quadrat ic 
l,x,x^ , (x~t)^ 




linear 
l,x, (x-t)^ 


Figure 6.1. A partial ordering of some spline spaces. 

Neither is a subspace of the other, however, so they cannot be formally 
compared (via testing). A solution of a sort is available if the 
quadratic and the linear are compared, since, as can be seen from the 

figure, the cP linear basis generates a subspace of the C® quadratic 
s pace . 

In general, a test to compare spline orders can be made bet%«en splines 

k-3 

of order k and k-1, both having continuity C . In the case of cubic 
splines, for example, we could allow continuity of the function and its 
first (but not second) derivative in order to determine whether the order 
could be reduced from 4 to 3 or increased from 3 to 4. Since a C'’ cubic 
has sufficient smoothness (at least to the eye), the procedure is not so 
objectionable. Considerably less satisfactory, however, are the cases for 
linear and quadratic splines. In comparing splines of order 3 and 2 as seen 
in Figure 6.1, the quadratic spline would be continuous but not its first 
derivative while in comparing splines of order 1 and 2, the linear spline 
would not even be continuous. Of course, the results of formal tests can be 
used in combination with informal comparison between SSE's of the models of 
interest to decide upon an acceptable model, and we recommend this approach. 

A backward elimination FORTRAN program using B-splines has been devel- 
oped for the purpose of reducing spline order using the nesting of some 
"sub-optimal" spaces as described above. Details for the appropriate B- 
spline hypothesis tests are given in reference 12. The listing, documen- 
tation and flowchart for the program are given in the Appendix, and we 
illustrate its use with the Indy data. While some statistical software 
packages could undoubtedly be used by defining "+" functions as in knot 




selection, no attempt was made to use them in this context. However, using 
the results of Section 4 as a guide, we surmise that several backward elimi- 
nation procedures would be suspect %ihile most forward selection algorithms 
should give fairly accurate results. Again, tolerance levels may have to be 
made small in order to force entry of certain terms. 

Figure 6.2 gives the FORTRAN program output for stepdown order selec- 
tion on the Indy data starting with a cubic spline (order 4) with the two 
knots in mid-WI and WWII as in Section 3. The progran compares splines of 
different order with the same continuity conditions, though other fits are 
given for information purposes. For this example, order reduction is made 
from cubic to quadratic to linear. Estimates of the B-spline coefficients 
and their standard errors are given for the spline of lowest order which can 
adequately fit the data, and the highest continuity conditions are imposed. 
For this case it is the Cr linear. 

A graphical display of these results is quite helpful, and Figure 6.3 
shows a partial ordering of the relevant spline spaces along with hypothesis 
test results and SSE's from the program. The dotted lines indicate the 
stepdown comparisons we wish to make, while the solid lines indicate those 
we can actually make through formal comparisons (tests). The importance of 
user input into the variable selection process is becoming more widely 
recognized, and here especially, because the formal tests available are not 
exactly what we would like. Consequently, we recommend the use not only of 
the formal tests, but also of informal comparisons between SSE's (or MSE's) 
of competing models using a display such as Figure 6.3. 

We illustrate this technique by going through Figure 6.3 step by step, 
and we will discover some interesting characteristics of splines along the 
way. We first observe that while a formal test is not possible between the 
C? cubic and the quadratic, it would not even be nec-'ssary since the 

cf quadratic has a smaller SSE than the cubic. A better fit is thus 

obtained with a lower degree! This phenomenon could never happen with poly- 
nomials, but such are the vagaries of splines. An informal comparison in 
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OF POOR 


TVE 3100T>€ST SPLIfC OF ORDCR K- 4 WITM MRXIttJI CONTINUITY C 2 
HPS SSE- 3B5.2SeeeZ18 AND tCC* 7.B63862S0 

CPN ORDER K- 4 KITH SUB-MRXIrtJM CONTINUITY C 1 

BE REDUCED TO ORDER K« 3 UITH mxinn CONTINUTIY Cl? 


YES. 

FOR K- 4 AND C 1 

FfABLE VALUE • 4.00 8 00000 OBSERVED F> 1.27322583 

SSE* 340.34966200 NSE* 7.41169494 


TVC SMOOTVCST SPLI« OF ORDER K- 3 WITVl MAXIMJ1 CONTINUITY C 1 
HRS SSE* 376.65994630 AND tCE- 7.53319093 

CAN ORDER K* 3 UITH SUB-MAXIMH CONTIHJITY C 0 

BE REDUCED TO ORDER K* 2 WITH rnxlMJM CONTINUTIY C 0 ? 


YES. 

FOR K* 3 AND C 0 

FTABLE VALUE * 4.00000000 OBSERVED F* 3.69172563 

SSE* 360.45371260 «SE* 7.67611901 


1>€ SnOOTTCST SPLirC OF ORDER K* 2 WITH imxiftjl CrNTINUITY C 0 
HAS SSE* 453.460O0B46 AND tCE* 8.89153115 

CAN ORDER K* 2 WITH SUB-NAXirtJI CONTINUITY C-1 

BE REDUCED TO ORDER K* 1 UITH MAXINLH CONTINUTIY C-1 ? 


NO. 

FOR K* 2 AND C-1 

FTABLE VALLE * 4.00000000 OBSERVED F* 285.46925596 

SSE* 328.52409313 NSE* 6.70457333 


PROCEDLRE TERMINATES. 


PROCEDLRE TEWINATES UITH L* 3; K- 2; C 0 


N COEF 

1 73.61675883 

2 00.49577629 

3 114.97961700 

4 157.47046102 


ST. ERR. 
2.19888373 
1.11724005 
.94655257 
1.10901397 


I***—******— FURTVER ITFORHRTION *»***»»» 
FOR K* 1 AND C-1 

SSE* 6070. 372 77 25 2 rSE* 116.7379379 3 


Figure 6.2. Output for order reductiion. Indy data. 
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385 


376 


453 


6070 


Figure 6.3. Partial ordering of spline spaces including SSE's 
and results of order reduction tests. Indy data. 

in going from the quadratic to the C® linear reveals an increase of 

77 in the SSE. While this increase cannot be formally judged insignificant, 
we may wish tc draw such a conclusion based on the results of the formal 
test which compares the C° quadratic with the linear: the larger 

increase of 85 is insignificant in that case. Having thus "safely" arrived 
at the C® linear, we must decide whether to further lower the order. The 
very large F value (285) from the program output which compares the C“^ 
linear and the C"^ constant splines reveals the importance of the linear 
trend. The big increase of 5742 in SSE from the C~^ linear to the 
constant is thus highly significant, and since the increase of 5617 from the 
C° linear to the constant (the desired comparison) is only slightly 

smaller, we conclude that the use of a constant spline fit is inadvisable. 

7. SPLINES IN SEVERAL VARIABLES 

A mathematical theory for splines in several variables is still devel- 
oping, and a "satisfactory" basis even in two variables has not been found. 
However, tensor products of either "+" functions or B-splines can be used to 
form a spline basis in several variables. While a tensor product basis is 
somewhat clumsy and its interpretation difficult, we explain here some theo- 
retical aspects of its use for the two variable case and give an example. 
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As in the example in Section 5, a spline-modified Taylor series expan- 
sion can be used to model aerodynamic force and moment coefficients. This 
time, however, Klein and Batterson (ref. 22) use splines in two variables, 
the angle of attack a and the sideslip angle b, to approximate the lat- 
eral force coefficient and the rolling and yawing moment coefficients. They 
use the yawing moment coefficient ag a typical example, and Cn 

can be expressed as 


C 

n 


C (a,b) 
n 0 

a 


6 

r 


+ C (a)p' + C r' 


p' - r' - 0 


+ C 


(a) 


6 + C (a) 6 

a n^ r 

r 


(7.1) 


where p and r are the rolling and yawing velocity and and 6r 

are the aileron and rudder deflection. They approximate the function Cjj 
(a,b) by 

ill „ 

C (a,b) - C + C b + I (A . + A,.b)(a - a.)'' 
n o 1 . , 01 li 1 + 

1*1 


t-2 
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B .(b 
oj 




tl 

I 

i*l 


l2 

Z D. .(b - b .) 

j-1 


(a - a . )'^ 
1 + 


(7.2) 


while the remaining functions in (7.1) are approximated by splines in a 
alone. Results from a stepwise regression using these terms are not as good 
as in the one-variable case, and some fine-tuning remains. 

From a theoretical point of view, the tensor product basis does not 
have the nice interpretation of knots and continuity constraints as in the 
one-variable case, even using functions. There is, however, a one-to- 
one correspondence between two-variable function terms and grid points, 
or nodes, and for this reason, we use the term node basis to refer to tensor 
products of the function basis. As before, we use right-continuous "+" 
functions so that 0° is 1. Tensor products of B-splines may also be used to 
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construct a basis for splines in two variables, and we shall see that the 

same advanta«*es and disadvantages of the one-variable case carry over. 

We discuss the simplest two-variable case in some detail: first order 

splines, i.e. step functions. Their application is somewhat limited, but 
there are several reasons for their detailed consideration. First and fore- 
most, splines in two variables are difficult to envision and manipulate, and 
consideration of the simplest case, namely constants, is thus highly desira- 
ble. Second, as seen in the example above and in Section 5, the estimation 
of aerodynamic force and moment coefficients using a spline-modified Taylor 
series expansion reveals the importance of using constants from both inter- 
pretative and numerical points of view. Finally, the two-dimensional cumu- 
lative distribution function is a first order spline in two variables. 

Thus, the constant case, while limited, has already shown its usefulness. 

We first discuss the node basis by way of example. Suppose breakpoints 
in the x variable occur at Xj, X 2 , X 3 and in the y variab’e at y^ 
and y 2 for data in Xg < x < x,^ and yg < y < y^. A function basis 
of order 1 in the x variable is (x - Xg)®,..., (x “ * 3 )^ and in the y 

variable is (y - Vg)®,..., (y " tensor product basis is formed 

by taking all the 4 x 3 * 12 products (x - x^)® (y - i “ 0 3; 

j ■ 0,..., 2. Each basis element in the two variables is thus a plane of 
height one bounded below by the line y * yj and on the left by the line x 
■ x^. Its support is thus a quadrant of a sort (l^l We call the 
intersection of these boundary lines, the corner of the quadrant, a node, 
denoted *ij. Figure 7.1 shows the relevant grid and nodes. Through any 



Nodes for a tensor product of functions. 


Figure 7.1. 




variable selection procedure, a model may be found whose terms are a subset 
of the 12 basis elements. Such a selection might result, for example, in 
the nodes shown in Figure 7.2 with the statistical model 

f(x,y) » Booxoyo^ + 302X0.^y2+ + 6lixi+yi + 

62 iX 2 ^yH. + B22*2+y2^. 632X3^y2^ ♦ e, 

where x. v. is an abbreviation for (x - x.)? (v - v.)+. We saw earlier 
the application of this technique to aerodynamic modeling. 



Figure 7.2. Nodes resulting after variable selection on 
a tensor product of functions. 

For splines of higher order, the same principles apply in forming the 
basis elements: they are the tensor product of one-variable functions. 

Knot multiplicities in ore variable result in node multiplicities in several 
variables. The absence or presence of a node or node multiplicity corres- 
ponds to the absence or presence of a certain basis element. There is thus 
some carry-over from the one-variable case in interpreting the role that 
basis elements play, and also in the fact that standard variable selection 
software may be used. The major drawback of this basis, as in the one- 
variable case, is computational. The basis elements do not have small 
support, so that roundoff errors get worse as computations increase. 

The computational difficulties present in the node basis lead to con- 
sideration of tensor product B-splines. While the formulation of the latter 
basis is straightforward, its interpretation and use in model selection 
through hypothesis tests are not. The polynomial degree and importance of 
knots in modeling are considerations that carry over from one to several 
variables, and unfortunately, so do their difficulties when using B-splines. 
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To compare differences in the two-variable case between the node basis 
and B-spline basis, we consider a simple grid with nodes indicated (*) in 
Figure 7.3. 


y?. 

yi 

yo 


if — *- 


Xq Xi X2 


Figure 7.3. Nodes for model (7.3) 


The statistical model for first order splines is thus 


f(x,y) = eooXQ+yo+ + Boixo^yi^ + Biixi^yi^ + e. 


(7.3) 


The function is a "true" spline in both variables except when yE[yg 
for then f is constant over [xg,X 2 ). If this model is represented with 
B-splines, each cell i is the support of a right-continuous plane which 
has height 1. Using the notation B^(x,y) for the basis element for each 
cell i, the model may be written 

4 

f(x,y) ® Y 6 .B.(x,y) + e subject to Bj * B 3 . 
i-1 ^ ^ 


This B-spline model is somewhat more complicated than the "+" function basis 
in its representation because of the model restrictions. It is also not 
obvious how to interpret the B-spline coefficients in terms of the presence 
or absence of nodes. 

These simple examples illustrate that the function terras are iden- 

tifiable and meaningful on a grid as nodes, just as they correspond to knots 
in the one-variable case. They thus hold an advantage over the tensor 
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product of B-splines from an interpretative point of view. As B-splines 
hold the computational edge, however, it would be desirable to identify the 
linear combinations of B-splines which correspond to tne presence or absence 
of nodes. The interpretation and use of tensor product splines of higher 
order is more difficult and remains to be examined in detail. 

8. SUGGESTIONS FOR FURTHER WORK 

There are potential research areas for both the univariate and multi- 
variate cases. In the univariate case, an efficient stepwise computer rou- 
tine using B-splines could be developed. This tiould give the user the 
choice of forward and backward procedures with a computationally efficient 
basis. The use of knot selection to fit data with loops could be investi- 
gated, and approaching the problem using the parsinetric technique of Smith, 
Price, and Howser (ref. 23), seems feasible. The successful use of splines 
in two variables has already been demonstrated (Section 7), but further work 
remains such as investigating fits to known underlying functions like we 
have done in the one-variable case. TVo-dimensional pictures in this case 
would be most helpful. Also, while the multivariate mathematical theory is 
still developing, interpretation of tensor-product bases from a statistical 
perspective could continue from that begun in Section 7. 
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APPENDIX 


Program Documentation 

Tvo FORTRAN programs have bt.en written which adapt stepdown procedures 
to B-spline regression. One program is for knot elimination trhile the other 
is for reducing the spline order. Theoretical details and appropriate ref- 
erences are given ' Sections 2 and 6. The programs are written in FORTRAN 
5 and have been implemented on both the ODU DEC-IO and the NASA/Langley CDC 
Cyber computers. Notation is patterned after that of de Boor (ref. 1), and 
definitions of parameters are given in the subroutine VL2NT, the second 
8ubrouti'’e called. All necessary input is read in or specified in subrou- 
ti'.e D.«” thf data, sample size NDATA, (initial) spline order K ■ degree 
+1, (initial) interior breakpoints and endpoints BREAK(*), number of conti- 
nuity conditions V(*) at the breakpoints, number of intervals L ■ # interior 
breakpoints +1, and tabled F value to be used in hypothesis tests. For 
equal spacing, the breakpoints and continuity conditions are most easily 
specified through a DO loo". Variables are dimensioned by one of three 
parameters (defined in comment staten^nts) which are specified in the PARA- 
METER statement at the beginning of the main program. 

Data must be interior to [bREAK(I), BREAK(L+1) ]. For tn ■ Indy data, 

X max - 61, so we arbitrarily set BREAK(l) * 0 and BREAK:l+ 1) - 62. V(I) is 
the number of continuity constraints at BREAK(I). For example, V(l) ■ 0 
means that the spline is discontinous at BREAK(l) while V(2) ~ 3 means there 
are 3 contiguous continuity conditions on the spline f at BR£AK(2), i.o., 
f, and f" are all continuous at BREAK(2). Note that V(I) must be 

less than or equal to K-l in order to have a "true" spline, not a polyno- 
mial, across BREAK(I). We always set V(l) ■ 0, though only for "symmetry" 
in the endpoint conditions, and V(L+l) need not be specified since it is 
never used nor referred to. 

The subroutine FLAG is designed to catch user input errors which would 
otherwise cause the program to terminate abnormally or give inacc.’.rate 
•■esults which may or may not be obvious to the user. Sample output detect- 
ing er ors in the input information of the Indy data is shown in Figure 
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OF POOR Q'JAl- . ^ 


TTC ORDER K • 4 

T»€ 4 WTERWfltS - • 3 
THE DPO6I0M N ■ 5 


BREn«>OINTS 





coNnmrrv coromore 

a 

4 

3 



IMDEX 

1 

2 

3 

4 


KEI«C 

D- 

4 

IGC( 

2)- 

4 


6 

7 

P 

9 


BREBKPOINTS MJST BE STRICTLY ITCREflSING. 

BREPKPOiNT as.sraaeeee is not i r<s; thw brerkpoint 
T>€ NUIBER of COKTINUITY conditions MUST be strictly 

LESS THRN !>€ SPH^€ ORDER K. V( 2)- 4 
AT BRElflKPOINT 33.50860000 IS TOO LRRGE, 

X VRUJE OUT OF RCNGE. 

X( L)> l.OWMMOOO IS NOT IN THE RPNGE EREPK(l)- 
TO BREPK(LRST)- 62.00000000 


7.SOB0000O 


5.00000 0 0 8 


STEPDOWi CfiftOT PROCEED. PfiOGRRM ABORTS. 


Figure A.l. Sample output detecting input errors. Indy data. 



S<*vor«l lint** in the programs are for plotting only. These are calls 
to the CDO system subroutines PSEUDO, INFOPLT, and CALPLT and the DO loop 10 
which calculates the spline values at the knots. 

For the knot elimination routine, input data and subsequently 
calculated information are printed by means of subroutines DATl and OUTNTS. 
This includes data values, spline order, number of intervals, dimension of 
the spline space, and knots. At each step of the procedure as indicated by 
the number of intervals L, the F-rat ios for the importance of each 
breakpoint are given along with the SSE and MSE. If a breakpoint can be 
eliminated, it is specified and the procedure continues to stepdown. If no 
breakpoint can be eliminated, the resulting number of intervals and spline 
order are given along with a list of the values of the B-spline coefficients 
and their standard errors. Sample ^utput appears in Figure 3.1, p. 8, in 
Section 3. 

As in the knot elimination program, the subroutines DATl and OUTNTS of 

the order reduction routine print input data and subsequently calculated 

information. In addition, at each step, the printout gives the SSE's and 

K“ 2 

MSE's for two splines of order K. one with continuity C and the other 
with continuity The hypothesis test is described in words with the 

results of the F test indicated. When further order reduction is not 
possible, estimates of the B-spline coefficients and their st.sndard errors 
are given for the spline of lowest acceptable order with highest continuity 
impv-ised. Additional information is given by including the SSE and MSE of 
the .lext lowest v>rder spline. Sample output for the Indy data appears in 
Figure 6.2, p. 41. 

Flowcharts are given in Figures A. 2 to A. 3 followed by the program 
listings. A full listing of the knot elimination program from a CDC is 
given, including the subroutines of de Boor (ref O that are used. For the 
order reduction progr.im we list only the main program and the subroutine 
SSHYP2, a variation of SSHYP appearing in the first program. 
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Figure A. 3. (concluded). 








oooortnono 


KNOT ELIMINATION PROGRAM LISTING 


mxRpn )AjOT(irmr,(xrTnn',Tnae6-ouTPur,'mPE2B,T)RPE2i,TRPe22) 
STEPDCW4 FCR BREPKPOinr ELXMINAnon FOR PIXO) ORDER K. 

7>c FijNcnoN pro ns first k-^ OERivATives rtJST BE corinruous. 

ronpK IS lerst t>c spmx sz2e, ^OR'm. 

mm IS AT LEAST N. UXTH MPKITUI COTfriTiJITY COTOITICMS, 

N-L>*K-1. urm NO OOMIH1I7Y COfOITIOre, N-Lm. 
tCITfOK IS AT least KaN. 


ApRArcTERcmRK-ieB.rfMAK-aee.KTmpK-aBeB) 

REAL BOOCF(mPK), Q<K7mP»(), DIAC(KTmP»0, TCrOMPX) 

* ,D03EF(mm), BRnmRK), BLf(mm),F(rerwo 
« ,CTRP6T(mPDO,AA(mA»<,mPD<), eRR0R(Nm<) 

« , tCE. N5H, SE(mA»<).FRAnO(mAD<) 

• .BB<mR»<,mAM), Lirf/(mA»(,mR»o,FB(mA»<) 

INTEGCR ERRDF,Mr,V,KQO(mP »0 

comm roATA. xcrooK), vcrcnm), ftable 
( xmm ^Apmw/ SREAKcmpK), cocrctarmo, l, k. vcrmo 

looim-e 

C er<TER DATA 

CALL DATKlOOUfn 

C SET THE MOT SEQUDCE 

QALL ML2Nr(SREAK.L.K,V,T,N,KErO) 

C PRELlMINflRY OUTPUT 

CALL OUTTfTS(BRCAK,V.L,T,N,K,KEIO) 

OCCK ITWT DATA 

CALL FLA60FLAS.N) 

IFdFLAC .CO. 1) 00 TO 25 

CALL PSEUDO 

C TEST FOR COMTIMUOUS K-l-ST DERIVATIVE AT EAOl KNOT 
JDCRIV-K-l 
I FHIN-FTAar 

Lm-L-1 

C GET THE LEAST SQUAfCS FIT, I.E.. T>€ B-SPlItC COEFFICIENTS 
CALL LSTSQ.i(T.N,K,0,DIAC.BC0EF) 

C LSTSQl CALLS B5PLVB, BOfAC, AtO BOELV 

C GET SSE pro rSE 
ERRDFtOATA-N 

CALL BSPU>P(T,BCOEr,N.K.DIAC,BREAK,COCF,L) 

CALL ERR,2l(r,ERR0R,EARIlF,SSE,nSE) 



c 
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ERRL21 CPLLS PPVnjU kHICH CALLS INTiRV 

LP1H>1 
DO la 1*1, LPl 

PB( I) -PPVnjLK BREAK , COer. L, K, BR£AK( I ) , 0 ) 
10 COKHNUE 


■acnammaM PLOTS 1 1 i 

CALL irrOR.T<0,NDATfl,X.l.r,l,0.,62.,74.,158.,l.,9, 

* SHINDY DfiTfi.l.lHY,0,5..4.,.75,.7S) 

CALL ITCOPLT(0,ND«TR,X, 1,Y,1,0. ,62. ,74. , 150. , 1. ,9. 

* SHirCY DRTA,1,1HY,22.5.,4.,.75,.75) 

CALL irrOPLT(l.U11,BREAK(2),l,rB(2).l,0.,62.,74.,158.,l.,9, 

* 9HU€)Y DflTO.l.lHY,l,5.,4.,.75,.75) 


IF(L .rC. 1) GO TO 12 

kRrrE( 20 .ii) ssE.nsE 

11 rOfWAK//' SSE-',ri6.a.5X,'«SE-'.F16.8) 


GO TO 9 

C TEST IrtWrPMCE OF EACH BREAKPOINT 
12 MMTE(20.2) L,PTAaX,SSE,r6E 

2FORNAT(//^' L-M3,5X,'F-TffflLE VALLE IS',Fl6.0,/-/ 

* ' SSE-'.F16.8.5X,'MSE-',Fl6.a/'/' 

« * F-^TIOS ARE: BREAKPOINTS ARE') 

3 DO 5 II«2,L 

ID-II 

CALL CNTRSTCID. JDePIV,N,K,L,T,BREPK,KEND,BRT,BLF.XOEF 

* ,CTRAST) 

C CNTRST CALLS XONT 

CALL SSHYPtBCOEF.CTRAST, Q.K,N,BRT, VAR SSH,H&»,HEF) 

C SSHYP CALLS FORSUB 

FRATIO<II)««&i/NSE 
kRITE(20,4) FRATIOdD.BREAKdl) 

4 FQRnAT(2F16.a) 

IF(FRATIOdl) .GE. FMIN) GO TO 5 
FWIN-FRATIOdI) 

KNOT- 1 1 

5 CONTINUL 

IF (FT1IN .LT. ’fTABLE) GO .TO 7 
l«ITE(20,6) 

6 FORMAK/' NO BREAKPOINT CAN BE ELIMINATED') 

GO TO 9 

7 nilN-FRATIO(KNOT) 





OF F. 


V- 




l«ITE(2e,a) BREMC(ICiOT) 

a rOBIftTC//' BRE«<POINT'.r7.3, ' IS ELIMINPTO') 


C RELABEL KNOT SEQUENCE T P6 kCLL AS BREAK. KEND. 
CALL REKN0T(KEND,KN0T,N,K,L,T,V,BREAK) 


AND V 


GO TO 1 

C PRINT RESULTING COEFICIENTS, STANDARD DEVIATIONS. AND F-VALUES 
9 CALL STDERRCQ.BCOEF.K.N.L.riSE.DIAG.AA.SE.LINV) 

C STDERR CALLS BCHINV AND MATVEC. 


CALL CALPLT(0.,0..999) 


25 STOP 
END 

C INDY DATA 

SUBROUTIhC DATKICOCWT) 

COtnON STATDENTS /DATA/ AND /APPROX/ ARE USED. 

C 

C THIS SUBRCUTINE READS IN T>€ DATA AND GIVES THE NUMBER AND 
C PLACEMENT OF THE KNOTS FOR TVC FITTED SPLI^€. 

h a h a t c ier (rm<-i00. NDm<-200. ktnmax-2000) 

INTEGER V 
REPL V X 

COmON / DATA / NDATA. X(NDMAX). Y(NDMAX). FTABLE 
COmON / APPROX / BREAK(H«X). COEFCKTTmo, L. K 
« . V(rf1»0 

NDATA • 55 
^RITE(20.51 

5 FOf»WT(' INDY DATA’/.'* YEAR Y X') 

DO I I-l, NDATA 

READC21.4) YEAR.Y(I),X(I) 

4 F0RMAT(I4.1X,FV.3. 1X.F2.0) 
fcRITE(20,2) YEAR.YCD.Xd) 

2 F0»WAT(I4.1X,F7.3. 1X.F3.0) 

1 CONTINUE 

C GIVE 1>c ORDER K AND HUBER OF INTERVALS L 
K • 4 
L • 3 

FTABLE • a. 00 

C GIVE TVC BREAKPOINTS AND CONTIHIITY CONSTRAINTS 
BREAK(l) • 0. 

BREAK(2) • 7.5 
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BR£AK(3> - 33.5 
BR£PK(4) - 82. 

V(l) -0 
V(2) - 3 
V(3) ■ 3 

RETIFN 

END 

SUBROUTlrC VL3<r(BREPK,L,K,V,T,N,KEND) 

COmjTES T>€ KNOT SEQUENCE T RND DirCNSION N FROM 1>€ BREAKPOINT 
C SEQUENCE BREAK, GIVEN TV€ SPLI^€ ORDER K, TVE M/1BER OF INTER- 
C VALS L, AND TVC MJMBER OF CONTINUITY CONDITIONS V(I) AT BREAK 
C (I). 

C 

input <nmwn»)» 

C BREAK (1) BREAKCL+l) TIE BREAKPOINT SEQUENCE. 

C I T>C (EMBER OF INTERVALS. 

C K....TK: ORDER OF T>C SPLI(C. 

C V(2),...,V(U TVE NUMBER OF CONTINUITY CONSTRAINTS AT 

C BREAK(2),...,BREAK(L). 

C 

C Kn t mn l n n * OUTPUT «»«»»» 

C T(1),....T(N+K) TVC KNOT SEQUENCE. 

C N TVE DirENSION OF TVE SPLII^ SPACE OF ORDER K. 

C KEND(I)....TVE INDEX OF TVE LARGEST KNOT EQUAL TO BREAK(I) 

C 

Q m m mmm METHOD 

C T>€ FIRST K KNOTS ARE SET EQUAL TO BREAK(l). T>C KNOTS ARE 
C TVEN sequenced so that K - V(I) KNOTS ARE AT BREAK(I) WITH 
C KEND(I) EQUAL TO TVE INDEX OF TVE LARGEST KNOT AT BREAK(I). 

C N IS SET EQUAL TO KEND(U AND TVE LAST K KNOTS T(N+1), . . . 
c T(nho are set equal to BREAKCL+1). 

C 

INTEGER K,L,N, I,V(1).J,ISTART, ISTOP. KEND(l) 

REAL BREAK(l), T(l) 

C SET TVE FIRST K KNOTS EQUAL TO BREAK(l). 

DO I I • 1, K 

1 T(I) • BREAK(l) 

C 

C FIND TVC INDEX KEND(I) OF TVE LARGEST KNOT EQUAL tq BREAK(I). 
KEND(l) • K 
DO 2 I • 2. L 

2 KENDd) • KEND(I-l) K - V(I) 

C 

C SET T(KEND(I-1) + 1) -...■ T(KEND(D) • BREAK(I). 

DO 10 I • 2, L 

ISTART • KEND(I-l) +1 
ISTOP • KENDd) 

DO 11 J • ISTART, ISTOP 

11 T(J) - BREAK(I) 

10 continue 

N • KEND(L) 

c 

C SET TVE LAST K KNOTS EQUAL TO BREAK(L+1). 
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V '■ l' 


DO 20 I • 1, K 
20 T(N4>I) • SRCPKCL-t-l) 

Rrnjfw 

END 

SUBROUTirC ajTNTS(BREPK.V,L,T,N,K,KOO) 

THIS suBForrirc is tor cutputing oh.y. it outputs pljl 

CnjLlNC flRGLfCNTS « T T E R VLeNT HRS BEEh CRLLED. 

«««««« input and OUTPUT*»i««** 

K....TVC SPLirC ORDER 

I TVE NLT1SER OF INTERVRLS 

N T>€ DirENSION OF 1>C SPLI^€ SPRCE 

BREflK(l),...,BREPKrL+l) T>€ BREPKPOINT SEQUENCE 

V(1),...V(U T>€ NUMBER OF CONTINUITY CONSTRPINTS RT 

BREPK(1),...BRERK(L) 

T(1),...T(N)....T>€ KNOT SEQUENCE 

KENDd KEND(U INDEX OF TVE LRfiOEST KNOT EOURL TO 

BRERKCl BREPK(U) 

DIVENSION T(l), KETffid), BREflk(l) 

INTEGER V(l) 

URITE(20,40) K, L, N 

40 rORMRTC//' T>€ ORDER K - 13/'^' TVE ♦ INTERVRLS U 13, 

* //' T>€ DIPENSION N - 13) 

HJITE(20,41) 

41 FORMAT (//' BREAKPOINTS', T20, ' CONTINUITY CONDITIONS ’ ) 

DO 45 J • 1, L 

45 M9ITE(20,42) BREAK(J), V(J) 

42 rORMRT(Fl6.8,T30,I3) 

LPl • L + 1 

kRITE(20,43) BRERK(LPl) 

43 FCRMRTCFIS.S) 

4^ITE(20,8) 

a F0RPt»T(/-/’ T INDEX’) 

PfV • N + K 
’“OUNT • 1 
INTEX • 1 

14?ITE(20,5) T(l), INDEX 
5 rORMRT(F16.8. 5X. 13) 

DO 7 J • 2, PPK 

IF (T(J) .EQ. T(J-D) GO TO 50 

URITE(20,12) ICOUNT, KEND( ICOUNT) . T(J), J 

12 FORMRT(T30, 'KEND( ' . 13, ’ ) • ', I3/T1,F16.8,5X, 13) 

GO TO 13 

50 kRITE(20,9) T(J), J 

9 F0fiMRT(F16.8,5X,I3) 

GO TO 7 

13 ICOUNT • ICOUNT f 1 
7 CONTINJE 

RETURN 

end 

SUBROUTirC FURGdFLRG.N) 

C THIS SUBROUTirC CVECKS FOR 
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C <1) BREPKPOINTS UWCH fift£ NOT STRICTLY INCRERSINS; 

C (2) TOO MRNY CONTINUITY CONDITIONS; 
c y I eattmo twqn 9oi 

C (4) X V«JUE3 uUT OF RRNGE OF T>C FIRST fiND LOST BRERW>OINTS. 

C 

PRRRf€TD?( TfW- 100, TWW-200, IOmW-2000 ) 

INTEGER V 

COmON / ORTR / NDflTP, X(Nim<), YCNOmo, FTPBLE 
COrtlON / APPROX / B?EPK(rrWX),C0eF(KTN1PX),L,K,V(ht»O 

DO 1 I-1,L 
IPl-I+l 

IF(BR£PK(I) .(Z. BREPKCln)) GO TO 2 

1 CONTINUE 

GO TO 4 

2 bf!ITE(2e.3) BREAK (I), BREAK (IPl) 

3 FORWTC/' breakpoints ftJST BE STRICTLY INCREASING. V 

* ' BREAKP0INT',F16.8,»<, 'IS NOT LESS THAN BREAKPOINT’, 

« Fie. 8) 

IFLAG-1 

4 DO 5 I-1,L 

IF(V(I) .GE. K) GO TO 6 

5 CONTIHJE 

GO TO 20 


6 l«ITE(2e,7) I,V(n,BREAK(I) 

7 F0RMAT(/' T>€ hi/IBER OF CONTINUITY CONDITIONS iHET BE STRICTLY'/ 

* 5X,' less THAN T>€ SPLIhC ORDER K. V( ' , 12, ' )•' , 12/ 

* 5X,' AT BREAKPOINT '.F16. 8, ' IS TOO LARGE. ' ) 

IFLAG-1 

20 IF (K .GT. 20? GO TO 8 
GO TO 10 

8 URITE(20,9) K 

9 FOfiMAT(/' K-'.I2, ' IS TOO LARGE.'/' TVE ORDER K MUST BE 20 OR', 

* ' LESS.’) 

IFLAG-1 

10 DO 11 I-l.NDATA 

IF(X(I) .LE. BREAK(I) .or. X(I) .GE. BREAK(L+1)> GO TO 12 

11 CONTirt£ 

GO TO 14 

12 I4?ITE<20,13) I,X(I),BREAK(1),BREAK(L+1) 

13 rORMATC/' X VALUE OUT OF RANGE. '/• X(',l4, ’)-’,F16.8, 

« ’ IS NOT IN T>€ RANGE BREAKC 1 )• ' ,Fl6.8/ 

* 5X, 'TO BREAK(LAST).',F16.8) 



OF PC C 


. ; \ 


14 1F(N .CT. «»T«) GO TO 16 
GO TO 18 

16 MR1TE(20,1T) N, NMT« 

17 rO»WT</' T>C DITO610N N-M2, 

* ' IS GREATER THAN T>C SWLE SI2EM4, 

IFLAG-l 

18 KTN«t<>Hi 

irCNDWW .LT. NDRTA) GO TO 19 
GO TO 

19 WRITE(2B,21) NDflAX.NDPTA 

. FORMATC/' OCCK PflRRrCTER STATEJ>ENT. ’/5X, ' NDMAX-MS, 

* ' MUST NOT BE LESS THAN T>C MJMBER OF DATA POINTS', 

* 14,'.') 

IFU«-1 

22 IF (hr»< .LT. N) GO TO 23 


O TO 25 


23 kPITEi20,24) M 

24 FORMATC/' OCCK PARAMETER STATEmiT. '/SX, ' tWN MUST NOT BE', 

* ' LESS THAN N-',I5) 

IFLAG-1 

25 IF(KTNMAX .LT. KTN) GO TO 26 
GO TO 28 

26 URITE(20,27) KTN 

27 F0RMAT(/' CHECK PPRAfCTEP STATE^E^T. '/5X, ' KTHMAK ML^^T NOT’, 

* ' BE LESS THAN KTN-', 14 , ’. ') 
in.AG-1 


28 IFCIFEAG. EQ. 0) GO TO 30 
kRITE(2B,29) 


29 FORMAT (//'/’ 

* ' STEPDOtM CPNOT PROCEED. 


PROGRAM ABORTS, ') 


30 

END 


RETURN 
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SUBROUTI^€ LSTSOKT.N K.Q.DIftC.BCOEF) 

cflas BSPLVB. Boroc. bqslv 

c 

COrrOi STfiTOW OflTfl is used. 

c 

C TWIS IS fl MODIFICATION OF DC BOOR'S SUBROUTirC L2flPPR, 

C PAGE 255. n irWTS T.N.K, Firos TVE least squares 
C APPfMXIMATION TO TVC DATA USING M0R( ARRAYS Q AND DIAG. 

C AND OUTPUTS TVC B-SPLIVC COCFTICIENTS BCOCF. 

C 

PARArCTER(KMAX<2e. NDMAX*2ae) 

REAL BCOO^CN), DIAG(N), Q(K,N). T(l). BIATX(KfW<) 

CtmjN / DATA / NDATA. X(NDNAX), YChCMAX). FTABLE 
C 

DO 7 J-l.N 

BCOEF(J) • e. 

DO 7 I-l.K 
7 Qd.J) • e. 

left • K 

LEFTTK • 0 
DO 20 LL'l. NDATA 

C LOCATE LETT ST X(LL' IN (T(LF.iT) ,T(L£FTt-l ) ) 

10 IF (LEFT .EQ. N) GO TO 15 

IF (X(LU .LT. T(L£rr+l)) GO TO 15 

LEFT • LEFT>1 
LEFTW • LEFTTK-cl 

GO TO 10 

15 CALL 8SPLVB(T,K.l.X(LL),LErr,BIATX) 

00 20 MM-1,K 
Dki • BIATX(tn) 

J • LEFTM<+tt1 

BC0EF(J» • DM*Y<LL) ♦ BCOGF(J) 

I-l 

DO 20 JJ-M1.K 

* Q(I,J) • BIATX(JJ1«DW 0(1, 

20 I • 

CALL BOCAC(Q.K,N.D’AG> 

CALL BCHSLV(Q,K,N. BCOCF) 

:tet\jrn 

□ 0 

SUBROUTirC BSPLV8(T,JV(I04.INDEX.X,LErT.BIATX) 

C CALCULATES TVC VALUE OF ALL POSSIBLY NONZERO B-SPLIICS AT X OF ORDER 
C 

C JOUT«MAX(JHIOH, (J^1)«(INDCX-1)) 

C 

C WITH KNOT SEQUENCE T. 

C DC BOOR PAGE 134-135 


PARArCTER( JMAX-20) 

INTEGER INDEX. JNIGH. left, I.J.JPl 

REAL BIATX(JHIGH).T(1).X. DeLTAL(JMAX),DELTAR(JMAX), SAVED, T091 

C DllCTCKjN BIATX(JOUT), TiLEFT+JOUT) 

C CUTRENT FORTRAN STANDARD MAKES IT I^POoSIILE TO SPECIFY TVC LENGTN 
C OF T AND OF BIATX PRECISELY WITV«UT TVC INTRODUCTION OF OTVCMISE 


If 


< 





Otitv. 

OF PC-'- V.. ‘V 


C SUPERFLOUS OOOITIONPL ARGUTCMTS. 

OflTfi J/l/ 

C SflVE J.DELTPL.DCLTflR (VPLID IN FORTORN 77) 

C 

GOTO (ie.20), 

le J -1 

BIATX(1)>1. 

IF (J.GE.JHI04) GOTO 99 

C 

20 JPl-J+1 


DELT«?(J)-T(L£FT+J)-X 
DQ-TRL ( J ) *X-T ( LjEFT> 1- J ) 

SR^<CD* 0 . 

DO 26 1*1, J 

TD91-BIRTX( I )/(XLTRR( I )+DELTRL( JP1*I ) ) 
BIRTX( I ) -SRWaHDELTRRC I )*TE»1 
26 SflVO)-DeLTRL(JPl-I)*TERM 

BIRTX(JP1)-SM3) 


ircEX 


j-jpi 

IF (J.LT.JHIW) GOTO 20 

C 

99 RETURN 

Em 

suBROurirc bchfac (u, nbrnds, m«3u. dia>^) 
c 

C COrCTRUCTS T>€ CH0LE3CY FACTORIZATION C ■ L • D • L-TRRMSPOSE. 
C SEE DE BOOR P. 256 
C 


INTEGER ^eRrC6. rfiOU, I. INRX. J, JNRX. N 
RERL U( NBRNDS. r«0U), DIAGcrROU). RATIO 
IF ( r«OW .CT. 1 ) GO TO 9 

IF (Wa.l) .CT.0.) Wd.l) • l./14(l.l) 

I 


C STORE DIAGONAL OF C IN DU G. 

9 DO 10 N-l,r«0U 
10 DIAG(N) • H(l.N) 

C FACTORIZATION 

DO 20 N*l.hR0M 

IF(W(l,N)+DIAGiN) .GT. DIAG(N))GO TO 15 
DO 14 J<l,r«RND5 

14 U(J.N) « 0. 

GO TO 20 

15 U(1,N) • l./Wd.N) 

infix • MI^«( NBRNDS- l.r«CU - N) 

IF (infix .LT. 1 ) GO TO 20 

JNRX • INAX 
DO 18 I-l. infix 

RATIO • W(I+1,N)»W(1.N) 

DO 17 Jd.JnAX 

17 ua.rn-i) • wcj.rni) - ucj-d.N ) •ratio 
JNRX • infix - 1 

18 W(I*1,N) • RATIO 
70 CONTIMJE 



EK) 

SUBROJTIhC BCH5LV(U,r«PNDS,r«0U,B) 

C SOLVES T>€ LirCPR SYSTOI C*X-B OF ORDER t«0W FOR X 
C PROVIDED U CONTflIhS T>€ OCLESKY rfiCTORIZPTIOM FOR 7>C BPr«€3) (SY«- 
C rCTRIC) POSITIVE DEFINITE NfiTRIX C flS CONSTRUCTED IN TVC SUBROUTIME 
C BOhrRC (QUO VIDE). 

C DCBOOR PnZ 258 

INTEGER r«PND6.r«OU. J. JNPX.N.rCNDNl 
REPL u()«prcs,r«ow),B(r«ou) 

IF (rROW.GT.l) GOTO 21 

Ba)-B(i»w(i,i) 

RETl^RN 

C 

C FOWRRD SUBSTITUTION. SOLVE L*Y-B FOR Y. STORE IN B. 

21 r®«m-f«prcs-i 

DO 3B N«l.r«OU 

jmcNir«( rcrcm , r«ou-N ) 

IF (JNRK.LT.I) goto 30 

DO 25 J-l.JNRX 

25 B(J-*N?-B(J-Mi)-W(J-*-l,N)»B(N) 

30 CONTirtJE 

C 

C BPCKSUBSTITUTION. SOLVE L-TRRNSP.X-D**(-1)*Y FOR X, STORE IN B. 

DO 40 N-M?OU.l.-l 
B(N)«B(N)4W(1.N) 

JMPK«MI«( r»em , r«ow-N) 

IF (JNRX.LT.l) GOTO 40 

DO 35 J-l.JMRX 

35 B(N)-B(N)-«(J>l,N)»B(J-»N) 

40 CONTirtJE 

RETL^^ 

END 

SUBROUTl. BSRJP (T.BCOEF.N.K.SCRTCH. BREAK. COBF.U 
C GRL1.S BSPLVB 
C 

C CONVERTS ’"HE B-REPRESENTj^riCN T. BCOEF.N, K OF SOTE SPLirC INTO ITS 
C PP-REPRESENTATION BREAK, COEF, L. K. 

C DE BOOR PAGES 140-141 
PRRAfCTER(Kr«X*20) 

INTEGER K,L,N, I. J. JPI.KMJ.LEFT.LSOFAR 

REAL BCXF(N),BREAK(l),COeF(K.l),T(l), SCRTCH(k.K) 

«, BIATX(KMAX).DIFF,FKnj,Sm 

C OirCNSlON BREAK(L>15.COEF(K,L).T(r+tK) 

LSOF»TR»0 

BREAKa)"T(K) 

DO 50 UEFT»K,N 

C FIND TVE rCXT NONTRIVIAL KNOT irfTERVAL. 

if (Tcueft+d.eq.tujtt)) goto 50 

lsofar-lsofar+i 

BREAK(LSOFAR^l )-T(UErr*n 
IF (K.GT.D goto 9 

COEF ( 1 . LSOFAR )• BCOer < LEFT) 



GOTO 58 

STORE T>C K B-SPLirC COEFT'S RELEVANT TO CURRENT KNOT INTERVRL 
C IN SCRTOK . , X ) . 

S DO 18 I-1,K 

18 SCRTCH(I,l)-BCOEF(LEFr-K>l) 

C 

C FOR J-l K-1. COrPUTE 1>C K-J B-SPLirC COEFT'S RELEVflKT TO 

C OIWNT KNOT INTERVfiL FOR T>€ J-TN DCRIVPTIVE BY DIFFERENCING 

C those for T>€ (J-DST derivative, and store in SCRTCH(.,>1). 

DO 28 JP1-2.K 

KNJ-K-J 

FMU-FLOfiT(KMJ) 

DO 28 I-l.KNJ 

DIFF-T ( LEFT+I ) -T ( LEFT+I-KMJ ) 

IF (DIFF.GT.8) SCRTCHd, JP1)» 

• <(9CRTCH(I>l.J)-SCRTCHa,J))>'DIFF)»#T<nj 

28 continue 

C 

C FOR J-0 K-l, Fir® T>€ values AT T(UETT) QT TVC J>1 

C B- SPLirCS OF ORDER J+1 l*CSE 3LPP0RT CONTAINS TVC CURRENT 
C KNOT interval FROM THOSE OF ORDER J (IN BWOO. THEN CONBIhC 
C WITH TVE B-SFLirC COEFT'S (IN SCRTCM( . ,K-J) ) FOU® EPR.IER 
C TO CCmiTE T>€ (K-J-i:ST DERIVATIVE AT T(LEFT) OF T>€ GIVEN 
C SPLirC. 

C NOTE. IF T>C REPEATED CALLS TO BSPLVB ARE THOUGHT TO GEIERATE 
C TOO MCH OVEWCAD, T>CN REPLACE T»C FIRST CALL BY 
C BIATXd)-!. 

C AND T>€ SUBSEQUENT CALL BY T>C STATEIENT 
C J-JPl-l 

C FOLLOkED BY A DIRECT COPY OF TVC LirCS 
C DEL’:5R(J)-T(LEFTd)-X 

C 

C BIATyOl) -SAVED 

C FPCN BSPLVB. DELTALCKmx' AND DELTAP'KnAX' WOULD HAVE TO 
C APPEAR IN A DlfCNSION STwTEmrr. OF COURSE. 

C 

CALL BSn.VB(T. 1.1. TfLEFT). left. BIATX) 

COEF(K.LSOFAR) -SCRTCH( l.K) 

DO 30 JP1-2.K 

CALL BSPLVBCT. JPl,2.T(L£FT).LEFT.BIATX) 

KNJ-K->1-JP1 

9LN-0. 

DO 28 I-l.JPi 

29 SLN-BlATXa)»>SCRTCH(I,KMJ)^SUN 

30 CC€F(Knj.LSOFAR)-SUM 

50 CONTINJE 

L-LSOFAR 

return 
END 

SUB90UTIVC ERR-2HFTAU. ERROR, ERRDF, 


SSE, MSE) 



no on o no n n nnnnn 


OF Pv'OR Qv'ALnY 


C«XS SUBPRCWRRfI PPVfiLUaNTERV) 

THIS SUBROUTirC OmtTES T>C ERROfi SS «SD MS. IT IS P 
MOOiriED VERSION OF DE BOOR'S SUBROUTITC L2ERR, PPGE 261. 
MSE IS T>C OUTPUTED fCPN SQUPRED ERROR. 

ppRRrrrcR( rm<> 100, Nim<>2Qe. KTT«v«<>2Q0e ) 

INTEGER ERRDF, V 

REPL n-flua). ERROR(l), MSE. Y, X. BREPK, COEF 
OirCNSION rrPUcNDPTP), ERR0R(NDPTP^ 

COrrON / DPTP NDPTP, X(NDMPX), Y(NCMPX). FTPBLE 
COttlON / PPPROX /' BREPKirf«>0. COEFcKTmw), L. K 

• . v(rm<) 


SSE'O. 

DO 10 LL-1, NDPTP 

rrPL(LL) • PPVPLU(BREPK,COEF.L.K,X(LL^0) 

ERRCRiU.' • Y(LL) - ETii^jaX) 

10 SSE • SSE ♦ ERR0RclL)**2 
MSE • SSE/ERRCF 

RET\i?N 

END 

REPL RK:TI0N PPVPLU(BREPK,C0EF.L,K,X.JDERIV) 

CPLLS 'INTERV 

CPLOJLPTES VPLUE PT X OF JDERIV-TH DERIVPTIVE OF PP EOT FROM P»»-fiEPR 
INTEGER JDERIV.K.L. l.M.NDUTtTT 
REPL BREPKCL'.COETCK.D.X, FmjDR.H 
PPVPUJ-0. 
nWDR-K-JDERIV 

DERIVPTIVES OF ORDER K OR HIOCR PREPRE IDENTICPLLY ZERO. 
IF iFmJDR.UE.0) GOTO 99 

FIND INTO I OF LPRGEST BREAKPOINT TO T>C LEFT OF X. 

CPU. INTERV ( BREAK, L,X. I, NDtmr) 

EVtSLUPTE JDERIV-TH DERIVATIVE OF I-TH POLYNOMIPL PIECE PT X. 
H-X-BREAKtl' 

DO 10 M-K. JDERIv+l.-l 

PPVPLU- 1 PPVWLU FmJDR > ««+K:0Er ( M. 1 1 
10 FmjDR-rmjDR-i 

99 RETURN 

END 

SUBROUTINE INTERV ( XT, LXT,\, left, ^ rUPG) 

C COTRUTES L£FT.MPX(I.1.L£.I.LE.UXT.PND.XT(IKL£.X) 

C DE BOOR PPGE 

INTEGER LEFT, UXT.rFUPG, IHl, ILO, ISTEP,rilDIXi: 

REPL X,XTiLXT) 

DPTP ILO I'- 

C SAVE ILO vP valid FORTTJPN STPTEj'ENT in TVC rCW 197T STANDARD' 
IHI»IL0*1 
IF ilHI.LT.LXT' 

IF iX.GE.XTiLXT" 

IF (UXT.LE.ll 


GOTO 20 
GOTO 110 
GOTO 90 



ILO-LXT-I 

iHi-ua 

C 

2B IF (X.GE.XTdHD) GOTO 40 

IF (X.GE.xraLO)) GOTO 100 

C 

C w-NOW X.LT.XT(ILO). DECREASE ILO TO CPPTIFE X. 

ISTEP-1 

31 IHI-ILO 

ILO-IHI-ISTtP 

IF (ILO.LE.l) GOTO 35 

IF (X.GE.XT(ILO)) GOTO 50 

ISTTT»-ISTU>»2 

GOTO 31 

35 ILO-1 

IF (X.uT.XT(D) GOTO 90 

GOTO 50 

C X.GE.XTdHI). INCREPSE IHI TO CflPTURE X, 

40 ISTEP-1 

41 ILO- IHI 
IHI-ILCHISTIP 

IF (lHl.GE.ua) GOTO 45 

IF (X.LT.XTdHD) GOTO 50 

ISTEP-ISTEP«2 

GOTO 41 

45 IF (X.GE.XT(LXT)) GOTO 110 

IHl-LXT 
C 

C wwFOJ XT(ILO).LE. X.LT.XTdHD. NPRROU TVE INTtJJWPL. 

50 MIDOL£-dLO*IHD>^ 

IF ( MIDDLE. EQ. ILO) GOTO 100 

C note, it is PSSUrCD THPT MIDDLE-ILO in CflSE IHI-ILO+1. 

IF (X.LT.XT(MIDDLE)) goto 53 

ILO-MIDDLE 

GOTO 50 

53 IHI -MIDDLE 

GOTO 50 


Cw*»iSET OUTTMT fiMD RETIJ»». 


90 

lELPG— 1 
LEFT-1 

RETURN 

100 

rELPG-0 
LEFT- ILO 

RETURN 

110 

mpG-i 

LEFT-LXT 

RETURN 


END 

SUBROUTirC CKTRSTd, JDCRIV, N, K. L, T, BREflK, KEND, 

« BRT, BLF, DCOEF, CTRftST) 

CPLLS BCONT 
C 

C finds T>€ CONTRRST COOTICIENTS for TESTING CONTINUITY OF T>€ 



C JDERIV-TH DeRIVOTIVE OT 7>C SPLirC FtrCTION «T HREflK(I). 

C 

CiMMita* input 

C I rUIBER Of INTERVPLS 

C T(1)....,T(NN().. .T>€ KNOT SEQUEj-CE 

C I T>C INDEX OF f»C BREPKPOINT Of irfTEREST 

C BREflK(l),...flREPK(L+U TVC BREflKPOINT SEOENCE 

C JDERIV NOrfCGflTIVE INTEGER GIVING T>€ ORDER Of T>€ DERI- 

C VRTIVE TO BE EVPUUfiTED 

C KEND(1),...KEMD(L) INDE5< Of T>€ U3RGEST KNOT EQUfiL TO 

C BREPK(l) BREPK(U 

C N. . . .DI^ENSI0N Of SPLI^€ SPflCE 
C K.... ORDER Of SPLII« 

C BRT, as, DCOEF...WORK flRRRYS Of LENGTH N 
C 

C ««»*« » OUTPUT 

C CTRRSTd CTRPST(N) THE CONTRfiST COEFTICIENTS USED TO 

C TEST CONTINUITY Of T>€ JDERIV-TH 

C DERIVfiTIVE PT BREPK(I) 

C 

CwoKUaeti method 

C TVC FTJHCTION SUBPROGRPtl BCONT IS USED TO COmjTE TVC VPLUE Of 
C TVC LEFT PND RIGHT LIMITS Of TVC JDERIV-TH DERIVPTIVE Of 
C RELEVfiNT B-SPLirCS PT aREPK(I). 

C 

INTEGER KEND(l) 

BRT(l). BLF(l). CTRPSTd). T(l), 

« BREPK(l), DCOEf(l) 

DO 20 JJ • I, N 
20 KOEF(JJ) • 0. 


DO 10 J • 1. N 
DCOET(J) • 1. 


COrfVTE VPLUE rOR RIGHT CONTINUITY 

If (KENDd'-K+l .LE. J .PND. J .LE.KENDd)) GO TO 30 
BRT(J) - 0. 

GO TO 40 

30 BRT(J) • BCONT(T,KOEF.N,K,BREPK(I),KENDd), 

* JDERIV) 


COMPUTE VPLUE fOR LEFT CONTINUITY 

40 IF(KENDd-l)-K+l .LE. J .PND. J .LE. 
BLF(J) - 0. 


KEND(I-D) GO TO 50 
GO TO 60 


50 BLF(J) • BCONT(T,DCO£F,N.K,BREPKd), 

* KENDd-1), JDERIV) 

C 

COMPUTE DIFFERENCE Of TVC LEFT PND RIGHT VPLUES 
60 CTRPST(J) ■ BRT(J) - BLF(J) 

DCOEF(J) • 0. 

10 CONTINUE 
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REPL RfCriOh BCONT(T,BCOCF,N,K,X,I,JtCRIV) 

CflLCUL«TlS VfiLUE PT X OF JTERIV-TH DERIVPTI'.t; OF SPLirC FROM B-REP. 

■mis IS P MODIFIED VERSION OF DE BOOR'S SUBROUTINE BVPLUE, 

PPGE 144. T>€ 0N.Y DIFTERENCE IS THPT THE LEFT-HRND KNOT 
IhOEX I IS irRUTED RPTVCR THPN FOUND IN INTERV. CONSE- 
QUENTLY. LII« 10 IS MODIFIED TO ir*>UT I PND LirCS 710 PND 
720 PRE OMITTED. THE ^LRPOSE IS TO PLLOW EVPLUPTION PT 
BREPKPOINTS WITH LEFT (0R RIOfT) CONTINUITY. 

PPRPrCTERC KJ1PX-20 ) 

INTEGER JDERIV.K.N, I. ILO, IM<, J. JC, JCMIN. JCMPX, JJ.KMJ.KNl.rTLPG 

* .rni 

REPL BCOEF(l),T(l),X, PJ<kMP9<),DL(KMPX),DRtKMPX),FKNJ 
DirCNSION TCN+K) 

BCONT-0. 

IF (JDERIV.QE.K) GOTO 99 

IF K-l (PT® JDERIV-0), BCONT-BCOEF( I ) . 

KMl-K-1 

IF (W11.GT.0) GOTO 1 

BCONT-BCOEF(I) 

GOTO 99 

«■ STORE THE K B-SPLINE COEFFICIENTS RELEVENT FOR THE KNOT INTERVPL 
(TCD.TCI+l)) IN PJ(l),...,PJ(K) PND COmjTE DUD'X-Td+l-J). 

0R(J)-T(14.J)-X, J-1 K-l. SET PW OF THE PJ NOT OBTPIHPBJE 

FROM imjT TO ZERO. SET PNY T.S NOT OBTPINPaX EQUPL TO T<1) OR 
TO T(NH0 PPPROPRIPTELY. 

1 JCMIN-1 
IfK-I-K 

IF (IMC.GE.0) GOTO 8 

JOlIN-l-IfK 
DO 5 J-1. I 

5 DL(J)-X-T(I+1-J) 

DO 6 J-I.KMl 

PJ(K-J)-0. 

6 DL(J)-DL(I) 

GOTO 10 

8 DO 9 J-l.KNl 

9 DL(J)*X-T(I+1-J) 

10 JCMPX-K 

rtil-N-I 

IF (rtU.GE.0) GOTO 18 

JCMPX-KHfll 
DO 15 J-l.JCMPX 

15 DR(J)-T(I+J)-X 
DO 16 J-JOWX.KNl 

PJ(J'fl)-0. 

16 DR(J)-DR(JCNPN) 

GOTO 20 

18 DO 19 J-l.KMl 


1 



19 DR(J)-T(I+J)-X 
C 

20 DO 21 JC-JOlIN,JCm< 

21 ftJ(JC)-BCCO-(HK+JC) 

C 

C **» OIFFOTNCE THE COEFTICIEKTS JDERIV TII1ES. 

ir (JDCRIV.EQ.e) GOTO 30 

DO 23 J-l.JDERIV 
Knj«K-J 

Fwn.OftT(Knj) 

ILO-KHJ 

DO 23 JJ-l.KWJ 

RJ( JJ)-( (ftj( JJ+1)-AJ( JJ) )/(DL(ILO)+DR( JJ) ) )»fWU 
23 ILO-lLO-1 

C 

C OmjTE VfiLUE AT X IN (T(I),T(I>D) OF JDERIV-TH DG?IVATIVE, 

C GIVEN ITS RELEVENT B-SPLirC COEFS IN AJ( 1 ) , . . . , AJ(K-JD£3?IV) . 

30 IF (JDERIV.EQ.KNl) GOTO 39 

DO 33 J-JDERIV>l,t011 

KMJ-K-J '■ 

ILO-KMJ 

DO 33 JJ-l.KMJ 

AJ C J J ) • ( AJ ( J J-*- 1 ) •DC ( ILO ) ( J J ) *DR ( J J ) ) / ( DL ( lUO ) ♦DR < J J ) ) 


33 

ILO-ILO-1 


39 

BCCNT-AJ(l) 


99 

END 

RETIRN 


C THIS IS FOR 1 DF WPOTVCSES. 

SUBROUTirC SSHyP(BC0E7.CTRAST,W,NBArffiS,N.PVAR,VAR, 

• SSH, MSH. HDF) 

CALLS rORSUB 

Q 

C finds T>€ variance of a C0NTT?AST and TVE re for TESTING THAT 
C TV€ CONTRAST IS ZERO. 

C 

c***«**imput»*»*** 

C LINV 1>€ INVERSE OF L OBTAItED FROM B C H I N V 

C CTRAST TV€ contrast vector OBTAirCD FROM C N T R 5 T 

C BCOEF....TVE B-SPLINE COEmCIENTS 

C W T>€ MATRIX FROM B C H F A C CONTAINING D-INVERSE 

C NBAtSDS EQUALS K 

C N TVC rUMBER OF ELEMENTS IN TVC CONTRAST VECTOR— 


C also T>€ DireNSION OF T>€ SPLirC SPACE 

C AVAR U0»< VECTOR OF LENGTH N EQUAL TO T>€ F=RODUCT 

C W(1,.)«A, I.E. D-INV*LINV*CTRAST 

C 

C***«**OUTPUT»****« 

C VAR. . . .TVC COEFFICIENT OF SIGMA-SQUARED IN TVC VARIANCE OF TVC 
C CONTRAST, I.E, TVC PRODUCT CTRAST-TRANSPOSE*LlNV- 

C TRArePOSE*D-lNV*t.INV»CTRAST 

C SSH, MSH, HDF.... TVC SS, MS, AND DF FOR TVC HYPOTVCSIS 
C 



ooooonoooooooonn 


or/r ' 
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C******MeTHOD*****’* 

C TVC PRODUCT LINV»CTRAST IS OBTAIhCD TVOt PREMJLTIPLIED BY 
C D-INV TVCN THRT RESULT IS PREMJLTIPLIED BY (LINV»CTRPST)- 

c transpose 

C 

INTEGER HDT 
REAL NUM, MSH 

REAL CTRAST(l), U(NBBNDS,N), PVAR(N) 

REAL BCOETd) 

Hjn- 0. 

DO 3 II-l.N 

3 HJ1 - MJ1 (CTRAST(II)*BCOeF(II)) 

CALL rORSUB(U,CTRAST,NBANDS,N) 

DO 1 1I*1,N 

1 PVAR(II) • W(1,II)*CTRAST(II) 

'VAR • 0. 

DO 2 II-l.N 

2 VAR • VAR + CTRAST(II)*PVAR(II) 

SSH • iNUri»»2)/VAR 
MSH ■ SSH 
HDF • 1 

RETLFN 
END 

SUBROUTirC FORSUE(U,AA,NBRNDS,r*?OU) 

SOLVES LY-AA FOR Y AND STORES IN AA 


»*»*»*iihput****** 

U. . .A MATRIX FED IN FROM B C H F A C AND CONTAINING IN ITS ROWS 
T>€ diagonal? or a P. D. SYTttTRIC MATRIX C 
NBANDS. . .TV€ BANDWIDTH OF C 
MWW. . .1>€ ORD OF C 

AA. . .TVE VECTOR OF LENGTH ?<?CW CONTAINING THE RIGHT HAND SIDE 


**»***0UTPIJT'****«* 

AA. . .T>€ SECTOR OF LENGTH M?OW CONTAINING THE SOLUTION 


*****«nETHOD**’***’* 

T>C FORWARD SUBSTITUTION ROUTirt FROM DCBOOR'S B06LV IS USED 

REAL MCNBANDS, rROW). AA(MTOW) 

IF (MJOW.GT.l) GO TO 21 
AA(1)-AA(1)«W(1.1) 

RETURN 

21 NBNDMl -NBANDS- 1 

DO 30 N-l,r*50W 

JMAX-riIN0(NflNDMl , MWW-N) 

IF iJMAX.LT.l) GO TO 30 
DO 25 J-l.MAX 

25 AA( J+N)-AA( J>N)-g( J-fl,N)*AA<N) 

30 CONTIHJE 
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Rmj94 

END 

SUBROUTirC REKNOT(KEISD,KNOT,N,K,UT,V.BREPK) 
RELABELS TVC KNOT SEQUENCE T( ) BY OmTTINC T>C LEAST 
SIO^IFICANT KNOT, BREAK (KNOT) 


KEND(I)...TVC INDEX OF 1>C LAE?GEST KNOT EQUAL TO BREAK(I) 
KNOT. . . INDEX OF T>€ BREAKPOINT TO BE OtIITTO 
N. ..DirtTSION OF TVC (OLD) SPLIIC SPACE 
K... ORDER OF TVC SPLIhC 
T. . . KNOT SEQUENCE 

V(I)..NUNBER OF CONTIHJITY CONDITIONS AT BREPK(I) 

BREAK... BREAKPOINT SEQUENCE 


N. ..DIttNSION OF (fCW) SPLirC SPACE WITH BREPK(KNOT) ONITTED 
Ta)..T(N)...(hCW) KNOT SEQUENCE WITH BREAK(KNOT) OMITTED 


SINCE BREPK(KNOT)'T(KEND(KNOT-l)+l)-...-T(KEND(KNOT)), kC 
RELABEL ALL T'S BEYOND. 


DI^CNSI0N KEND(l), T(l), BRE.CK'1) 
INTE(3ER V(l) 

Il-KEND(KNOT-l)-H 

I2-KEI®(KN0T)-t-l 

Jl»NH<-I2-t-l 
DO 1 KT-l.Jl 
Kl-KT-l 

1 T(I1+K1)«T(I2+K1) 
N»N-(K-V(KNOT)) 

DO 2 II-KNOT.L 
BREAK( II) -BREAK (II-«-l) 

IF(II .EQ. U 00 TO 2 

V(II)-V(II+l) 

KEND ( 1 1 ) -KEND ( 1 1- 1 ) -HC-V ( 1 1 ) 

2 CONTIMJE 

L-L-1 


RETURN 

end 

SUBROUTIhC STDEPR(W,BC0EF.K,N.L.MSE.BB.AA.SE,LINV) 

ILLS BCHINY AND NATVEC 

THIS SUBROUTirC COTPUTES TVC STANDARD ERRORS OF TVC B-SPLItC 
COEFTICIENTS AND OUTPLfTS noi. 

REAL W(K,N),BCOEr(N),MSE,BB(N,N).SE(N).Llt^(N.N),AA(N,N) 
CALL BCHINV(W,K,N,Lir^O 
VRITE(2B,10) L,K 

10 rORMAT(///’ PROCEDURE TERMINATES WITH L-M3, ' AND K-M3// 



N 


COEF 



* ' 


S.E. ') 


DO 11 11*1, N 
DO 11 J>1,N 

li BB(JJ,II)-wa,II)»<-INV<II,JJ) 

CflLL MRTVEC(N,N,N.BB,LINV,AP) 


DO 13 H-l.N 

SCai)*SQRT(PP(II, 

H?ITE(20,12) II,BCOeF(II),SE(II) 

12 rORMRT(13,2F16.8) 

13 C0NTIH£ 

RETURN 

END 

SUBRelUTl^€ BCHINV (W, NBfiNDS, f«0W, INV) 

C finds L-INVERSE l*€RE L IS THE LOiCR TRIflNGLLflR MRTRIX 
C IN T>€ CHOLESKY FftCTORIZfiTION OT 1>C BPNDED SYtICTRIC P.D. 

C MPTRIX C flS CONSTRUCTED IN TVC SUBROUTirC B C H F fi C. 

C SEE DC BOOR. P. 256 
C 

C«aM[4Dmi INPUT o xinaa * 

C r«OW IS TIC ORDER OF TIC MPTRIX C. 

C NBRNDS IS TIC BPNDWIDTH OF C. 

C W CONTRINS TIC CHOLESKY FflCTORIZPTION OF C PS OUTPUT 

C FROM SUBROUTIIC B C H F P C WITH ROtfi 2 THROUGH NBPNDS-1 
C CONTPININS THE NON-ZERO PND NON-UNIT DIPGONPL ENTRIES 

C OF L. 

C 

Onom OUTPUT 

C INV TIC INVERSE OF L. 

C 

C METHOD ***** * 

C TIC LIICPR SYSTEM L*L-INVERSE - IDENTITY IS SOL'CD FOR 
C L*INVERSE BY SUCCESSIVELY FINDING THE COLUTTC OF L-INVERSE 
C USING TIC FOajPRD SUBSTITUTION ROUTIIC IN B C H S L V. 

C 

INTEGER NBPNDS, HROW. J, JMPX.N.NBNDNl 
REPL U( NBPNDS, HROU), INV(HROU,r<ROU) 

IF ( HWW .GT. 1 ) GO TO 21 

INV(l.l) • 1. 

RETURN 

C 

C STORE TIC IDENTITY MPTRIX IN INV 
21 DO IB J-l.HWW 

DO 10 I-l,H?OW 

IF (I .EQ. J) GO TO 20 

INV(I.J) • 0. 

GO TO 10 

20 INV(I.J) • 1. 

10 CONTINUE 

C 

C N4J USE FORUPRD SUBSTITUTION FROM B C H S L V. 
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25 

30 

40 


NENDMl > NBCr<DS - 1 
DO 40 J-1.M90U 
DO 30 N-l,r«CM 

JMBX • WN0(NBNDm,N?OW-N) 
IF (JMPK .LT. l) 

DO 25 I • l.JMRX 

INV(14«,J) - INV(I+«,J) 
CONTINUE 
CONTIHJE 


GO TO 30 

W<I+1,N)«INV(N,J) 


O® 

SUBROUTINE MfiTVEC(N,l^,M,X.Y,2) 
FINDS 1>€ rWI MfiTRIX OR VECTOR Z l+tlCH 
X»Y UM35E X IS rWN AND Y IS r<W1. 


RETURN 


IS THE F=RODUCT 


REAL X(N.I^), Y(M1,M), Z(N.M) 

DO 1 I • l.N 
DO £ J • 1,N 
Z(I,J) - 0. 

DO 3 K-l,h#1 

3 Z(I,J) • 2(1, J) + X(I,K)*Y(K, J) 

2 CONTIK£ 

1 CONTINUE 


RETLW( 


nooooooono 


cr?-;. 


i- 




ORDER REDUCTION PROGRAM LISTING 


PROGRPT »=^T^(I^FUT. OUTPUT, TflPE6«OUTPUT,TPPE20.TflPE21' 
STEPDOI.W fOR REDUCING SPLIfC ORDER <rOR «X INTERVALS 
SirULTRhCOUSLY) i^ILE KEEPING 1>C KNOTS FIXED PND PGSUniNG 
K-2 CONTINUITY CONDITIO^S 


NEMRX IS PT least T>€ SPfPLE SIZE, NDATA. 

rm< IS PT least n. with mpkinjh continuity cotctraints, 
N-L+K-1. UI7VI NO continuity CONSTRAINTS, N-L*K. 
KTNIPX IS PT LEAST K*N. 

PPfiPrCTER (mPX-130,NDMPX-200,t<TmW*2000) 

REAL BC0EF(N1PX),Q(KTN1P>O,DlAG(KTm»O,T(Nrm<) 

* ,LINV(KTTm<),OCOCF(rr«><),BRT(h*1PX),aE(N1PX) 

* ,PP(NMPX).VPR(MW<),B<Nt;><),C(N1A>0,PTRP(rmo 

* .r(NDMPK),E3?ROR(NDm'<),MSH.nSE,£Z(NnPO<) 

* ,KMPT(MiPK,rNp>o,rB(rriPX) 

* ,wvpR(rtttx),cc(NiPX),a ^^m<),cTRAST(^m<) 

INTEGER ERRtF,HDF,V,KEND(N1AX) 

COrtlON /DATA/ NDPTP,X(NDf1AX),Y(NimO,FTABLE 
COrWON /APPROX/ BREAK(fWX),C0Er(KTT«PX),L,K,V(TI1AX) 

ICOUNT.0 
C ENTER DATA 

CALL DPTKICOINT) 

C GET TVC knot sequence 

CALL VL2NT(a?EAK,L,K,V,T,N,KEND) 

C PRELIrUNPRY OUTPUT 

CALL OUmrS(BREAK,V,L.T,N,K,KEND) 

OECK irPUT DPTP 
IFLPG-d 

CALL FLAG(IFLAG.N) 

IFdFLAG .EQ. 1) GO TO 25 

CALL PSEUDO 

IEND-0 

Lm«L-i 


C LC UILL TEST -HAT T>€ K-l-ST DERIVATIVE IS ZERO IN ALL INTERVALS 
C 

C GET THE least SQUARES FIT. I.E., TVC B-SPLINE COEFFICIENTS. 

1 CALL LSTSQKT.N.K.Q.DIAG.BCOEF) 

C LSTSQl CALLS BSPLVB, BOCAC, AND BCHSLV. 
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C GET SSE flHD rSE 

ERRDF-NDflTf»-N 

CPU. BSPLPP(T,BCO£F.N,K.DIftG.BREPK,COer.U 
CPU. ERRL21(F,ERR0R.ERRDf.SSE.MSE) 

C ERRL21 CPLLS PfMtU WICH CPU.S INTERV. 

iraEKD .EQ. 2) GO TO ae 

LPl-L-t-1 
DO 10 I-l.LPl 

10 FB(I)-PPVPLU(BREPK,COeF.L.k.BREPKa),0) 

C i t i l n t, rfT l n l n M 0 > nnn« n.0T5 trtnnnr^irnnnrti 

CPU. irfOPLT(0.NDPTA,X, 1.'-. 1,0. ,62. .74. . 158. , 1. , 

* 9.5HINDY DflTA.l,lHY,0,5..4.,.75,.75) 

CPU. lhFOR.T(0.NWTP,X. 1. Y. 1.0. ,62. .74. , 150. . 1. . 

* 9. SHINDY DATA.l IHY, 22.5.. 4... 75.. 75) 

CPU. IrrOR.T(l.U11,BREPK(2).l.FB(2).l,0.,62.,74..158..1., 

* 3.9HINDY DATA,1,1HY.1.5.,4.,.75. .75) 


W11-K-1 

KM2-K-2 

KM3-K-3 

irclEND EP.^' GO TO 8 
IF(V(2).EC Xj-C) GO TO 12 
WRITE(20,r) K.KTE.SSE.MSE 
15 FORNPTC//' 

» //' tVE SMOOTHEST SPLirC OF ORDER K-M2.2X, 

* 'WITH MAXIHJM COf.TIHJITY CM2,2X/^, 

» 'HAS SSr»',F16.8,2X, 'PND nSE-',F16.8) 

IF(K.EQ 1) Gi> TO 8 
RITE(20,16) X,KM3,Km,Krt3 

1 . F0RNPT(/' CPN ORDER K- ' . 12, 2X. 'WITH SUB-MAXIMLTi CONTINUl'iY C', 
I2.2X/SX, 'BE REDUCED TO ORDER K-'. 12, 2X. 

* 'UITH nR^.*tri CONTIHJTIY C',I2. ’ ’' ) 

DO 18 II-2.L 

le V(II)*K-2 

CPU VL2NT(BREPK,L ’'.V.T.N.KEND) 

GO TO 1 


C TEST FOR LCXC3 ORDER WITH THE HYP01>€SIS MATRIX KMAT. 

12 DO 90 JJ*1,N 
80 DCOEF(JJ) • 0. 

DO 2 III-l.N 
DCOEFdll) • 

DO 61 Il-l.L 

KMPT ( ! 1 , 1 1 1 ) • 9CONT ( T . DCOO-' , M , K . BREAK ( 1 1 ) , KEND ( A I ) . KM 1 ) 
M-(II-1)HS^I^ 

31 CT!M)-KMPT(IT.,III) 

DCOEFCTII) . 0. 
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2 CONTIHJE 

CRLL SSHVP2(BC0EF.CT,Q.K.L.N,flA.DIA6.VfiR,SSH.Ma<.«r.B.C 

• ,«TRP,UVflR,CC,CnWST) 

C SSHVP2 CflLLS rORSUB (M) nftrvEC. 

nwTio-fSi/teE 

IF(FWTIO.GE.rr«a£) go to 5 
kMITE(20.3) 

3 FCRViUy/’ ns. •) 

W?ITE(20,31) K,KM3.FT«BL£.nWTIO,SSE,MSE 

31 FOfWPTC FOR K-M2.2X, 'AND CM2.2X/5X, 

• 'FTAaX value • '.F16. a. 5X. 'OBSERVED F-',F16.9 

• 'SSE-'.F16.a.2X. 'MSE-'.F16.8) 

K-K-1 

CALL VL^^■(BREAK,L.K,V,T,N,KE^C) 

GO TO 1 

5 IE»-1 
«ITE(20,6) 

6 F0»«T(x/* NO. ') 

t«ITEcaB.31) K.KM3.FT«BL£.nWTIO,S3E.nSE 
M7ITE(2B.32) 

32 F0*WAT(//* PROCEDURE TEWIINRTES • •nrm o n ota 
DO 71 II-2.L 

71 V(II)-K-1 

CALL VL2<T(BREAK.L,K.V.T.N.MEJ®) 

GO TO 1 

C PRINT REaJLTING COEFFICIENTS AND STANDARD ERRORS. 

8 U?ITE(20,13) L,K,KM2 
1? '^ru/y PROCEDURE TERMINATES WITW l-M2. CM2 
//' N COEF ST. ERR.') 

CRLL STDERR(Q,BCOEF,K,N.L.MSE.DIAG.AA.SE,LINV) 

C S TDERR CAIXS BCHINV AND MPTVEC. 

IF(K .EQ. 1) GO TO 25 

IEND-2 

K-K-l 

DO 40 II«2,L 
40 V(II)-K-. 

CALL V4_2NT(BREAK,L,K,V,T.N.KEND) 

GO TO I 

30 U?ITE(20,29) KMl,KN3.SSE,riSE 

29 FORIPTC/'/'' «*»»»«»»»»»»»»»» RJ?T>CR irfORMATION *»*»»*«'/ 

• 5X, 'FOR K-’.I2,2X, 'AND C’.IZ.ZX.' 

• 5X. 'SSE»',F16.8,2<, 'MSE-'.F16.8) 

CALL CALH-T(0.,0.,999) 

25 STOP 
END 
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OF FOC;. 


SUBRa^■I^€ SSHYP2(BC0e=-,CT.W.NBfM)S.NC0ri.N.A 
• .P''(W.VfiR,S^^,^BH,HDF,B,C,ATRP,^MW,CC.CTRfiST) 

INTEGER HCT, NCON,N,NBRNDS 
REPL S9< 

REPL PVW(NCON.N).fl(NCON,N),U(NBflHDS,N),CT(l) 

REflL VPR(NCON.NCON), B(N), C(NCON), PTRP(N.NCCN) 

REPL MSH, UVfiR(NCON,NCON). CTRPSTCNCON.N), BCOET(N) ,CC(N) 
DO 1 I-l.NCON 
DO 2 JJ-l.N 

M»JJ+(I-1)*N 

CC(JJ)-CT(M) 

2 CTRRST(I,JJ)-CC(JJ) 

CPLL rORSUB(U.CC.NBPHDS.N) 

DO 3 J-l.N 

3 ft(I,J)-CC(J) 

1 CONTIHJE 

DO 4 rr-i,NC0N 
DO 4 JJ-l.N 

4 PVPR(II,JJ)-Wa,JI)*P<II,JJ) 

DO 5 I«l,N 

DO 5 J-l.NCON 

5 ftTRP(I.J)-P(J,i: 

CPLL MRTVECCNCON.N.NCON.PVPR.PTRP.VPR) 

DO 6 I-l.NCON 

m-NCOtt-I+l 
DO 6 j-i.m 

6 WVPR(I,J)-VPR(I+J-1.J) 

CPLL BOFPC (UVPR.NCas.NCON.DIPG) 

CPLL nPTVECCNCON.N. l,i.TRPST.BCOeF.B) 

CPLL rORSUB(UVPR.B,NCON,NCON) 

DO 8 J-l.NCON 
a C(J)»UVPR(1.J)*B(J) 

CPU NPTVECd.NCON.l.B.C.SSH) 

NSH-SSH^ON 

HDT-NCON 

RETURN 

END 
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